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ABSTRACT . 

Seismic refraction techniques were employed in crustal 
Studies of South-Central Saskatchewan.A unique arrangement 
of the sources and receivers was set up to obtain three 
dimensional structure over a 300 km equilateral triangle.The 
interpretation of the data presented here is based on the 
three shot-reversed sides of the triangle. | 

The evaluation of the collected data was undertaken by 
use of several techniques including ray tracing,inversion of 
head waves and synthetic seismogram simulation.The derived 
models reveal a particularly complex many-layered structure 
of the crust below the study area.Layers are gently dipping 
the dip angle never being greater than 2.5°.Low velocity 
zones-probably associated with high COnduCctiwaity 
material-have been detected below all three profiles. P-wave 
velocity in these layers is close to 6 km/s.Significant 
lateral heterogeneity is found below the southeasterly and 
northeasterly profiles (B and C).There is good evidence for 
artavirequstawest of Pine «A. “The fault. strike ts’ ©S.6&. 
passing through the earthquake prone area of Bengough.Upper 
mantle velocity is practically invariable under the study 
area and close to an average value of 8.13 km/s.Depth of the 
Mohorovicirc discontinuity varies from 36 to.45 km.oThe crust 


thins eastward and thickens southward. 
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1. INTRODUCTION 

In August ISS be the COL CRUST sgroup ‘conducted =1 es sthird 
cooperative seismic survey in south-central Saskatchewan.The 
1981 program consisted of four survey lines A,B,C and D .The 
bocatron of the™protiles issshownein™ figure 14. “Refraction 
and wide-angle reflection techniques were employed.There 
were two main objectives for the experiment.The first aim 
was to obtain information about the structure of the earth's 
crust underlying the northcentral segment of the 
intracratonic Williston basin.In particular line A_ was 
placed close to the NACP conductivity zone in hope of 
obtaining some seismic expression of this anomalous area. 
Secondly the experiment tested a novel type of crustal 
refraction data gathering method.In this system,following a 
suggestion of E.R.Kanasewich, three profiles(lines A,B and 
C)were arranged in an equilateral triangular form so that 
three dimensional information was obtained by the 
combination of inline reversed shots and also by the 
broadside recordings from opposite corners of the 
triangle.This thesis deals with the interpretation of the 
three in-line recordings of the triangle.The analysis of the 
collected seismic data was undertaken by the use of three 
different methods: 1)Inversion of first arriving head wave 
“intormation,2)Ray “tracing "of wide angle reflected energy, 
and 3) The interpretation of amplitude information using 
synthetic seismograms.Future analysis of the overlapping 


Mroacside' data ima the central portion-of the triangle jwill 


ai) 


E nee | ind cu 
a, a bre Wee 
. 
soitedstee - 

oor ee 


e'isaas eit to gurtoedge eds, tyoma” nobtanys® aks 


ais te  sneneee-. teirasteion eH 4 piiytas 
5 atii “sledisweg aT plese. oogeriery vines ” itl 
to sepa ong aos Gavi sSaoS ROA ate az ‘es hats 


molpgrons ehaAgy Pd nsteesuens  sinaiee once 


onieollot moteyve eidze hi vbote $a bo ivedane 8285 nod ao 
hfe 8A saint Dlasbitew eerie , 48) wae LUA » AS 10 pape Ou 


A+ Ae meee % af ine weP Qe tex re HH ti ‘iia né ne SONS7Ts 


pe 
G 
r 


bencedde wav rbisveisotal | fenotensmd a een 


eiz wi oal@e SOAs s30%e “Sptvevyat entini - noisenle sailed 
. Bein : , ase ees. 
beconge moi? spnibsenes | ENSS 


aes, 


S 

~ 

4 

a 
“a 


eas he oly a: teaquaint of? daw Alesh etveds Bier. ignel * se 


git to ateylone od? elpnetay ont Ye apnibierss oni ot 
er . p 


geniz.?o ee. wit ye gedergaeiae Say £288 senten! ap 
i f 7 


os 


sveuw beat odtedate tages) To ie tevewnt (t seat 


ae a 


vb tee nated Mvecanagee shaw Bey ee aa 
oviey oocsmagetal . ene shot iigma to anaes 


t ; i 
oniqggalteve ef9- te elaybans old 
v- a Vs 


iliw elpapiae saa 3e). Noh s 208, isianen 
| a Oia 


, 
oo2 vys 


J > p> er > 
SS Ky perm pare 


AzOT 


i 


5 > >>>> > >> >> >> >>> > >>>> 


APOT x VE k 
vs VA ee 3 
ve leet olensel silat alot! lela, Selalalnalalaleal- LE 
fi TT] fe d x ae ce 
Pa . r Lan fee ee 
es er pwn, <n wt 
"e ee a - 
“9 a oar 


be 
M90T + ie 
ene ee ue 


ASOT 


/ 
B 
gee 
\ 
Oo xeon, \ 
ori rs, 
\ 


“ 
Ahaunaron 


m Great Valls 


HMETT 


=Fi igure 1.... Profiles x B ismic 
ate pos : POC UROL Ene 1981 CO-GRUST seismic 
© 
cS 
wv 


47N 


er 
] 
ae | 
— a as 


e ~~ : ue 
i» | Le 
4 ‘a 7 
1 5 a ’ vi a hd 
‘ a a 
; a 
| a | 
Hin > = 
oes 
/ eRe a 
| | = a’ 
| | ae F 
lad : of oc’ e* 
) rs 
| i x i 
| | ; 
| r a 
4 
| . . ’ 
. OU Bo cs A 
| ' : ‘ 
1 ¢ ~ 
} ' “ ® 
; i id ; 
iu . een 
i a a a . 
i had ' 


, 
a 


ee ras 


See ae ‘ xSP 


: oy 


i” a y 
iat j 


be based on the models presented here. 

The inferred crustal models are compared with previous 
seismic studies in eastern and western Saskatchewan . The 
Operatucnadiegend theoretical Saspects Pot. the pworke are 


discussed. 


b = ey mando stand 
rs oye wes 


t 


ee De i ae von 


“ g.: 
hi 


ng Bae 
Me 


Set 


2. GEOLOGY AND GEOPHYSICS OF THE STUDY AREA 


2.1 Geology of the study area 

Geologically Canada is divided into a number of 
regions, (figure 2). The area of the 1981 experiment is the 
central-south Saskatchewan. It is placed entirely within the 
Interior Platform which is the part of the North American 
craton that subsided slightly in the past and now is covered 
with Phanerozoic sedimentary rocks. To the north and east 
the craton crops out forming the Precambrian Canadian shield 
which 1s the largest continuously exposed area of 
Precambrian rocks in the world. To the west of the Interior 
Platform the craton is bordered by the Cordilleran orogen. 
The southern boundary of the western Canada sedimentary 
basin can be defined only in an arbitrary manner. If all of 
the Williston basin is included, the boundary is probably 
best drawn by the Sioux uplift in southwestern Minnesota to 
the Black Hills and thence to the Little Rocky mountains and 
westward to the Lewis re 

By far the greater part of the basin is floored by 
rocks which were intruded or metamorphosed 1800 million 
years ‘agowe The @span-cof-1200vm.y4, prior™to Cambrian time 
permitted deep erosion and the development of a marked 
degree of tectonic stability within the interior plains 
area. 

Contours drawn on the surface of the Precambrian 


indicate the existence of three basins and two arches in the 
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Interior Plains area, between latitudes 49 and 60°. From 
northwest to southeast the basement features have been 
named: 1)Northern Alberta basin, 2)Peace River arch, 
3)central Alberta basin,4)Sweetgrass arch and 5)Williston 
basin. It is in this latter area that the 1981 experiment 
was entirely placed. The landscape of the Williston basin is 
Simple. Near surface sediments of Cretaceous age are 
essentially flat lying and have an average elevation of a 
few hundred meters above sea level.The older and gently 
Gem nl aiescepeawias. beens spurred by) deposits#of glacial 
till, sand and gravel. The Phanerozoic sedimentary section 
tS? Wp. Gon 3.0 km thick)’ near vy Southern ‘Saskatehwean.. The 
sediments thin northward to the exposed craton of northern 
Saskatchewan, (figure 3). 

This difference is due to varied basin subsidence and 
sedimentation in different regions.Subsidence was 
accentuated at various times throughout the 
Phanerozoic.Basin sedimentation was locally affected by 
movements of structural regions or blocks of the Precambrian 
basement delimited by lineaments as well as by solution of 
massive halite beds which also resulted in the collapse of 
overlying strata. The base of the stratigraphic column is a 
succesion of early Paleozoic clastics followed by Ordovician 
to Mississipian carbonates and evaporites and Jurassic to 
Pleistocene clastics.The clastic rocks of the Cambrian lie 


unconformably on the basement, figure 4. 
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EFTgureys.... basement configuration in western Canada (from 
Geologic history of Canada,1965). 
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Figure 4.... Reference geological column for southern 
Saskatchewan indicating order of geological systems.Arrows 
ingitate Major Unconformities. . 
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Geologic information on the Precambrian comes from 
studies of the exposed shield (Bell,1971) and studies from 
core samples(Burwash and Culbert,1976). The basement under 
the Interior Plains appears to be a south-western 
continuation of the exposed Precambrian shield. In 
particular, the composition of the basement under. the 
Williston basin shows a high concentration of metamorphosed 
and volcanic and metavolcanic rocks.Percentage of 


Sedimentary rocks 1S about 35%. 


2.2 Previous seismic studies 

Previous crustal seismic studies in Saskatchewan 
include the works of Chandra and Cumming(1972) and Green and 
co-workers (1980). Their interpretation deals with the 
Suffield-Swift Current profvte end the | Weyburn—-Brandon 
profile respectively. Their results are illustrated in 
figure Se.) Both works indicatel: a relatively thick 
multilayered crust at the margin of the Williston Basin. The 
Riel discontinuity (Clowes,Kanasewich and Cumming,1968B) of 
tesent in both models. The Weyburn Brandon profile has been 
reinterpreted by Moon and De Landro (1981),figure 6. Their 
model does not show evidence of the Riel discontinuity. 
Instead they introduced a positive velocity gradient in the 
lower crust. There is also the Leross Oungre profile (B) of 
Pacem CO> CROs) mes eexper ment. wiuch has beensincrerpre ced) by 
Kazmierczak, (1980), figure 7.Finally the (unreversed) line 


C of the 1979 survey was interpreted by M.Shahriar 
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Crustal models from refraction experiments in 


Figure 5 


rea 


cn Albe 


he 


Southern Saskatchewan and Sout 


Manitoba, 


H 


7 2 
> ? 1 of z= [ ) 
% ° - =4 . 4 — —— ee 
- - * 
; ~ wa = 
* ' saad ee - 
Se i weet | 7 
’ ‘ ? 1? ae 
t aad i 


eL.aetsler? Suadsgiads 


> = 
F . 
| |e 
Hl =. oae* —, 
LP ere mis b 4d 
at =» rs we 
: | ee 
t O08 
t * 
hal 
| ~ - ~ an 
- <ibd 
ae = 
J ~ 
c —- — 
4 vi 4 
ed 
> Soke = - Vv a - as 


am 
a oP 


él 00% miss ROGS 287 os: mew elshoa fubousd 


esteailkA Arehitua’ San PANES RIES GlSe 


PO eG 
% a 
¥ 


B AB 


| 
CHURCHILL 


I 
RIOR 
7 ISUPERI 
PROVINCE : | PROVINCE 
1 ct 
km/sec ¢- boundary -? 
zone 
ks le en BL | 
Zee ee «O97 
TOO > Sie) 
13.00 ee (- 0.007) 
6.49 | 6.50 s/sec 
\ 
grad (0.007) l grad. (0.010) 
| 
| 41.00 
45.00 i ee ek: 8 1 a 


Prequeei ies. . Crustal model from CO-CRUST 1977 east-west 
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(1982).The derived horizontal model is illustrated in table 
Herheeproriie runs east-west ana latitude’ -of 749 °)" 938". 
from Limerick (longitude 106°, 44') to Carlyle (longitude 
102°, sDeO VHence Mit crosses "our siine-A®tand ‘therefore = "rt 
will be a useful reference for comparison. The shot point 
was at the west end of the profile which was shot reversed 
(as line D) during the 1981 survey. 

Furthermore ,combined refraction and reflection 
interpretations revealed the existence of a major crustal 
Paure= =Witnine scene +€hurcihi Pl) province (Green et art 
1976,b,Kazmierczac,1980), at a longitude of 103 °. There is 
aevatiation rnrverustal ‘thickness*-across= the Siauliv® “of Wat 


least 5 km,the upthrown block being at its west side. 


Layer Velocity Thickness Depth 
km/sec km km 
1 Siz 2:8 0.0 
2 6.0 225 aed: 
3 632 630 5 20 
4 62.1 SiO dilpero 
5 6e5 520 14.0 
6 ae ti 620 on 
if 654 4.0 Zone 
8 Giri Hehe =: 2220 
9 ON = 44.5 
Table 1 
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3. DATA ACQUISITION AND PROCESSING 


3.1 Data acquisition and description of profiles 
Thennorthassouthwline, linetA,was elocated simethe tcentral 
eequon jofiithercrustal midcontinental conductivity “zone. This 
anomalous area was first described by Camfield and 
Gough(1977).The northeasterly line,line B,was placed between 
Wynyard and Swift Current.The south-easterly line,line C, 
extends from Swift Current to Regway at the Canadian-U.S. 
border.The total lengths of the lines were 287.9 km, 288.6 
km and 288.4 km respectively. Ideally,the experiment would 
cover the triangle ABC by only three shots. However one 
immediately realizes the di Daremlty concerning the 
availability of 180 individual stations' which would be 
required for such an experiment, let alone the huge cost of 
the operation. Only 45 stations were-available,adequate to 
cover only one profile at a time.Hence three shots were 
required for each line(two inline shots plus the broadside 
shot) ive. ra © totale of » 9 es shots siwere detonated for the 
triangle. All recording sites were selected along roads such 
that no location deviates from the staight line more than + 
ikm.The average separation between recording sites was 6.5 
km.Ideally the spacing would be around 1 Km (close to the 
wavelength of the seismic energy). Elevations,latitudes and 
longitudes of all recording sites and shot locations were 
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SOovdr nates isin the rangero:, 270.0005) r.e0 £5935, me sShot 
point to receiver distances were calculated using a geodetic 
computer program by supplying the geographical coordinates 
as input data.Time of shots,receiver coordinates and 
elevations as well as instrumentation and shot point 
information were compiled by Dr.Hajnal of the Department of 
Geological Sciences of the University of Saskatchewan. All 
recording systems were pre-programmed to a set of specific 
time windows. The digital interval of the data was 0.016675 
(60 Hz).Consequently the Nyqyist Esau ey was 30 Hz.The 
University of Alberta recorded with the "University of 
Toronto digital event recording system".Vertical as well as 
N-S and E-W horizontal component of motion were recorded.In 
addition,the absolute time signal(WWVB) was incorporated as 
a fourth channel in our magnetic tape data (see figure 
8).This was of particular importance in order to check the 
accuracy of the times of the first data samples for each 
Station.In addition to the above set of equipment we also 
Recorded Imtofa Hwanear array, Consisting of 12/7 T1=36;2 \Hz 
vertical component geophones with spacing of 268 m.Two 
traces of the array were included in the major line 
information. 

All 1981 refraction data were compiled in digital form 
on a 9-track, 1600 B.P.I.(bytes per inch) magnetic tape. 
Each block of data,-a total of 24 blocks per channel-in the 
Original master tape contained 5 seconds of data and 


occupied 4800 bytes of memory.The logical record length was 
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80 bytes, hence the format by which each block was written 
was 60*(5E16.6).The total set of the recorded data occupied 
436 files on the magnetic tape.Each file was preceded by its 
header information which consisted of ten care 
images,containing detailed information about the recording 
Station.We decided,for our data processing purposes,to 
reformat the data and we copied them into a new tape in 
which they were written in integer +*2 form, each block 
consisting eve 1200 data samples(hence a total of 6 
blocks/channel). The logical record length was 20 


bytes,hence the format of the new tape was 120*(10F8.0). 


3.2 Presentation of the data 

The first step towards the evalution of seismic 
refraction data is the display of the records.For this 
purpose,the form of record sections is generally chosen.A 
record section consists of seismic traces on a time-distance 
graph.The simplest way is to plot travel time versus 
distance wSing «shot instant “and location as the’ zero 
reference. Such a presentation requires a large page size or 
an inconveniently small time scale. The other form of 
display uses a reduced travel time instead of time. The 
reduced travel time t' is defined as t'=t-x/V where t is the 
fame from shot detonation,x is the offset distance and V is 
the reduction velocity. Such a presentation has the effect 
of rotating all the travel time branches clockwise around 


their point of intersection with the t axis,the amount of 
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rotation=~being @“anversely ‘proportional "*to'’ the ® reduction 
velocity V.The result is better time resolution over smaller 
Space.V was chosen to be equal to 6.5 km/sec for the data in 
CUS WOGKe et esl Giliicance? OL velco that "Mirco VeLOCiIty 15 
close to the apparent velocity of the waves of interest it 
will produce a display showing the head wave phases all 
arriving at the same reduced time.Hence this kind of display 
facilitates the correlation process. 

Since we planned to process only the P-wave phases, all 
OUr  recora Sections Show * only the “traces “for vertrcal 
component seismometers. It must be mentioned that no 
elevation corrections were applied to the data for two 
reasons. 1)Such corrections would require the knowledge of 
near surface velocities in detail.2)The amount of correction 
WOUlds De  mNSstqniticant (€.03 Sec maximum) since the 
differences between station elevations was small (less than 


Gur meters.) 


3.3 Power spectral analysis 

ASeC1sSmMre “trace 1s a. function of tame; f(t);representing 
the velocity, sof “the <gqround at the site of the sensor.Its 
Fourier transform f(w),yields the phase and amplitude 
Spectrum of the ground velocity.In many cases however it is 
more appropriate to consider the power spectrum of the 
Signal.The power spectral estimate is equal to the smoothed 
Square of the amplitude spectrum. 

The technique we applied for our power Spectral 


estimates was as follows.The periodogram was computed by use 
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of a discrete fast Fourier transform.An example is shown in 
figure 9.Clearly the periodogram needs to be smoothed and 
for this purpose we applied the Daniell (or box-car) window 
in the frequency domain according to the equation: 


m 
a 
2m+1] ) 


PW) = = F (W-J) P* (W-d) 


j=-m 
W=Ussl 2eerk. GN / 2. 

The Daniell window was chosen because 1)It gives non 
negative estimates and 2)It gives equal weight over the 
entire range of its frequency bandwidth. N'(the full length 
of the time series) was equal to 512,while N was 360 
(=6s).The time ees was padded with zeros for 
361<3<512.The parameter m which determines the resolution 
was taken equal to 3.The effect of the application of the 
Daniell window on the periodogram of figure 5 can be aes in 
figure 10, which shows our spectral estimate for trace 
(2471, Jomthne Stationiwasmat a distance of 145.5 km from the 
shot.in figure 11 and figure 12 one can see the power 
Snectralep estimacesmiror ijehe traces(10, 1,1) and) (36,1,1) ‘at 
offsets of 70.7km and 265.9 km respectively. By comparing 
the three spectra, one can see clearly the displacement of 
the seismic energy towards lower frequencies aS we move 
further away from the shot.This is expected of course if we 
recall that the amplitude A(x,t) of a one dimensional plane 
wave is given,apart from a geometrical spreading term,by the 
formila A(x,t)= Acexoli(k(w) + 1a(w))x-wt] wiere x is ‘the 


distance the seismic wave has travelled and a,the 
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attenuation coefficient, is proportional to the frequency of 
the wave.Thus, at large distances from the _ shot,high 
frequencies will be greatly attenuated. 

It is of interest at this point to see the power 
Spectrum characteristics of the background noise of the 
records.This will be useful later when bandpass filtering 
will be applied on the records.Thus, power spectra of the 
noise for several traces were computed;figure 13 and figure 
14 show two typical examples of such Spectral 
estimates.Careful comparison of the computed power spectra 
leads to the following conclusions: 

1% TestS on noise show that the energy of the background 
noise of the records is concentrated between 1.5 and 6 Hz. 
Zee In nearly all cases the energy of the noise is 
decreasing beyond 15 Hz. 

33 Spectral estimates of the signals themselves (after the 
onset of the seismic energy)- figures 10,11,12 show that the 
Significant energy is always between 1 and 17 Hz. 

The above conclusions suggest that a band of 5-17 Hz is the 
appropriate one when bandpass filtering is to be applied on 


the records. 


3.4 Filtering 

Broadly speaking, a filter is a device or a physical 
process that operates on a time history and changes the time 
history in some manner.The principal use of BYNCe rer nc git 


seismic data reduction is frequency smoothing of the data 


ew 


sy 


ae 
Yanevoes? 


$ 
Cr aed * a 
“4 4 
1 
- 
> 
a 
. 
. 
‘ 
jo “J 
‘ > 
\ 
” 
7? 


d 


mA 


fo 


\ 


~é 


= re hs || ant fe 
1 *, SUR. oe 


wong? : 


._¢ Thaw. 


= eval, wiz 


e98..9ms23 yeaten a a 


a 


ar 
4 
d a of { "2 
a ty FO £ ad 
% 
, 
apt i afi 
f ow 
he 
5 ) d.6 4784 22% 
— ~ < e 
« S 75.31 
¢ ¢ - - laa 
! — 


205.00 


00 


See Cee 
U1.00 82.00 heen OLE 


aude ee tle) Uv ss 


O00 


opaele 


oe 15.00 22.50 30. 98 
PREQUENCNS NIHZ 


TRACES nie ly TEST PON NGELSE 


Pigquee 3.05). 


Oise 


Speetnalmestamate ofetheinoise of record 


15) 


ee 
(eS) 
eae ee TES 15.00 22.50 
Panicl Ehlert seers 
EAC EaS Saleepeon CN NGLSE 


Pigucee| Se .cceteseeon noise fOr trace 6,377.1 


SOOO 


26 


Ss ~~ ae 
Si? 3243, 27000 ne Feet 
n . —- arin t 


27 


under analysis. This means that the filter is designed to 
pass signal frequencies in a desired band and to attenuate 
all other unwanted frequencies, outside the band. Digital 
recursive filters allow the» greatest flexibility —in 
frequency smoothing and in this study a recursive zero-phase 
bandpass’ Butterworth "filter "will be Puseful. “A detailed 
analysis of the eight-pole bandpass Butterworth filter can 
be found in Kanasewich (1981). A FORTRAN program for such a 
filter was written by D. Ganley in 1977 during his graduate 
Studies at the University of Alberta. As an example of the 
effectiveness of the filter,figure 15 compares a record of 
unfiltered data with the same record after filtering with 
two different bands, 5-17 Hz and 7-13 Hz . The improvement 
of the quality of the record is clearly visible.However care 
must be exercised in timing the onset of first energy with 
data filtered by a zero phase shift filter. Finally a band 
of 5-17 Hz was decided for our records following the 
conclusions of the power spectral analysis discussed in the 


previous section. 


3.5 Despiking the data 

As a general comment on the original data, the noise 
level was rather low, and the first breaks of the seismic 
energy quite clear and distinct from the noise background. 
However, a significant problem with the University for 
Toronto recorders was their small dynamic range leading to a 
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appearance was that of a series of spikes, whenever’ the 
amplitude of the signal was large. 

These rapid and high amplitude variations are not 
attributable to a physical process but rather to saturation 
effects, either saturation of the digitizer or saturation of 
thes lamplifiers® AASs ai) mactermor itacta the whighy sensitivity 
required to obtain good signals on early-arriving head waves 
made it very difficult to adjust the gains of the amplifiers 
during the experiment. | 

A computer program was written which attempts to obtain 
the correct sign of the signals while making no attempt at 
obtaining the correct amplitude. 

The program works as follows: The data are passed 
through several tests and accesses to these tests decided 
upon the operation to be applied, which generally consists 
Otay ce Ppropereavteration of che sign and theyampliuude sof the 
Samples. The program makes use of several empirical 
parameters, which were decided after careful examination of 
a large number of spikes and which concern fabd etn amplitude 
of a spike, differences in amplitudes between adjacent 
Samples etc inawivgure Wie. aren fshowns ie) “Top syot the 
page:record 32, shot 6,vertical and horizontal components 
(original data ), b) The same record after application of 
the despiking operation and c) The same despiked record 
Sitter aoolication oo: the Butterwortin filter band) 5-17  jH2-A 
floweichantweot sche despiking iprogramipis, shown in figure 17. 


Very little use was made of heavily spiked areas so _ the 
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interpretation does not depend critically on reliability of 


these data. 


3.6 Normalization 

The =coneept of normalization 1S quite  wmportant in 
order to obtain usable amplitude information from _ the 
seismic section. Different seismometer responses ,coupling 
between the seismometers and the ground and different gain 
Seu ngs rare atc ltmeeractors:) Ghat il may cause considerable 
amplitude differences from trace to trace. The purpose of 
normalization is to minimize the effect of these factors. 

It 1S expected theoretically that the average power of 
the records should decrease monotonically with 
distance.Thus,multiplication factors for each trace were 
calculated in order to make the power values of the traces 
decrease with distance ina linear fashion according to the 
equation POWER(J)=POWER(1)-B*DIST(J). J is the number of the 
Station, increasing with distance and DIST(J) its offset. 
The range within which the power values will vary,between 
the first and the last trace of the section, can be easily 
controlled by proper adjustment of the slope B.In figure 
18,(record section of line A,shot 1),maximum amplitude 
values are the same for all the traces,while figure 19 shows 
the same section after application of the above technique of 
normalization. The numbers at the end of the traces are the 
Salculated ss ampiitieation,.factors,..of the traces. It can be 


seen that a considerable improvement in relative amplitude 


information has been achieved. 
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4, SEISMOLOGY AND THE EARTH'S CRUST 


4.1 The main crustal wave groups 

In a layered elastic isotropic medium which is a first 
order valid approximation for the earth's crust,two kinds of 
body waves can exist:compressional P-waves and transverse,or 
shear S-waves.There are four main P-wave phases usually 
observed in seismic crustal studies:Pg, P*,Pm and Pn. 

a) The Pg phase is a wave penetrating into the upper 
part of the basement and is generally observable up to a 
distance of 60-100 km. It has a velocity between 5.5 and 6.4 
km/s.In a continuous profile it will comprise several 
prograde segments of travel time branches with the velocity 
MmcneaS kg ewich weCEStance I nhis sesplictinga natures otlthe Po 
group will not be detected if the spacing of the recording 
Stati Ons was stoou lance (=10..km..) sbhegquctands for “granitic” 
layer. 

b) The P* phase is a P-wave refracted through an 
intermediate layer in the crust with a velocity near 6.5 
km/s. The upper boundary of this layer has been called 
Conrad discontinuity in Europe.The P* phase is often the 
most prominent first arriving event in crustal seismic 
sections. 

c) The Pm is a wave reflected at the crust/mantle 
CoaAnS UC lOREeOG OldGamwy. ' tiers. vomnetrograde branch son 4a 
time-distance presentation which shows bee a strong 


velocity discontinuity exists between the crust and the 
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mantle.The Pm curve is generally well expressed by large 
amplitudes in the range of the critical point (where the 
angle of incidenceis i=arcsin(V(Pg)/V(Pm))). 

d). The Pn waves are diffracted from the top of the 
mantle. They have a linear relationship on a prograde branch 
on a graph of arrival time versus distance from the shot. At 
sufficient distances(between 150-300 km) the Pn is the first 
arrival wave,the amplitude being usually very small and its 
frequency Significantly lower than the Pm. Its velocity is 
between 7.8 and 8.3 km/s. The Pn phase has long been 
interpreted as a head wave.The observed amplitude,however,is 
usually greater than that predicted for head waves, implying 
that the velocity change at the MOHO is not exactly 
Step-like but it has a finite gradient within a transition 
zone.Currently,instead of interpreting the Pn as a head 
wave, a more likely explanation is given in terms of guided 
waves within a high Q layer several tens of kilometers in 
thickness at the top of the mantle. Anisotropy of the Pn 
velocity has been reported by several authors(Hess 


,1964,Rait et al.,1969) in an oceanic setting. 


4.2 The composition of the crust 

The great value of seismic velocity structures is that 
they may allow us to make an interpretation about the 
fectonwe history) ofssthe Treqion'@land*iott the ethology fat 
various eee Heat inks measurements of seismic velocities 
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pressures and temperature conditions simulating those in the 
crust can be compared with velocities derived from seismic 
criustalirecords tvominter? iniormat tons faboutie there kinds) i vot 
_ rocks present at various depths.Jeffrey's .(1937) early 
hypothesis that the crust varies from the top to the bottom 
between an acidic and intermediate composition is probably a 
Still valid approximation. However, modern work has'7 refined 
this oversimplified picture. 

Recent results from geological and geophysical studies 
can be summarized as follows. Crustal refraction 
Surveys encountering Precambrian crystalline basement rocks 
usually measure upper crust velocities from 5.8 to 6.1 km/s 
which are typical of light granitic rocks and metamorphosed 
sediments. The upper Srisit:, below the sedimentary 
cover,consists of layered metamorphic rocks such Se aoe sees 
and schists,often penetrated by granitic intrusions.The 
grade of metamorphism increases with depth. In the shallow 
or middle crustal refraction sections one also has slightly 
higher velocities which range from 6.1 to 6.4 km/s.These are 
Often ttattrabuted t to meqniatmeite rocks which are believed to 
comprise most of the middle crust and on a larger scale they 
may behave as intrusions into the overlying cover of 
metamorphic rocks(Haller,1958,Smithson et al. 1975). Seismic 
velocities from the lower: crust are most often in the range 
between 6.5-6.9 km/s.These velocities are too low for 
gabbroic rocks and there is strong evidence that they are 


associated with a lower crust predominantly composed of 
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granulite facies rocks (Smithson,Shive and Brown 
,1977).However,lower crustal velocities in the range of 
7-7.5 km/s have been reported in several seismic. studies 
showang rthatierparntsecof thesloweracrust fare truly gabbroic. 
Upper mantle velocities in the range of 7.8-8.4 km/s are 
certainly attributed to ultramafic rocks and most probably 
to dunite. 

Numerous attempts have been made recently(Burke and 
Dewey,1973, Gibb and Walcott ~197%+) to explain the 
development of the SCAR has through plate tectonics.The 


plate tectonics hypothesis implies a complex series of 


events affecting continents Eiciluding ml ines NG 
continent-continent collision and continent-arc 
collision.Horizontal movements of lithospheric plates 


predominate.Surely such processes must lead to distinctive 
geophysical sSignatures.The first most obvious feature that 
we would expect is fused continental blocks of different 
composition and therefore of different physical properties. 
In fact, the lateral variation or blocks that have been 
Suggested for continental crustal structure(Kosminskaya and 
Zverev,1968;Smithson and Brown,1977) may be, in part, a 
vestige of Precambrian plate tectonic processes.We would 
expect continental blocks of different composition to have 
different thicknesses since Precambrian areas should be 
approaching isostatic equilibrium.Discontinuities of great 
Verlabihityinandebstructural complexity s should be, located 


between these blocks. 
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Plate tectonics theory would explain the gently folded 
blocks generally identified in the lower crust as part of an 
early succession marked by horizontal tectonic movement 
whereas the more steeply dipping zones often detected in the 
upper crust are a later development during a lock-up stage 


promoted by collision. 


4.3 Low velocity zones 

Low velocity zones in the upper and middle crust have 
been postulated by Landisman and Mueller(1966) in north 
Germany,Pennsylvania ,Central Japan and east Montana,and in 
regions of the Canadian Shield by Braile(1977).According to 
Snell's law,no seismic rays originating from a near surface 
source can be returned to the surface from within a low 
velocity zone.Hence,an interpretation of the seismic data 
based on the first arrivals alone will give no evidence of 
the LVZ. Furthermore, one obtains depths greater than the 
true ones for all the lower boundaries. To recognize the 
existence ofarsuch” a’®izone:y one Sineeds'® "to identify othe 
reflections from the top and the bottom of the LVZ in a high 
quality crustal seismic section. The use of synthetic 
seismograms for this purpose is of great help. 

A Satisfactory explanation for the existence of 
velocity inversions in the crust has not yet been reached.It 
1s however an apparently common feature of the crust, its 
presence being observed in geologic sections that have 


undergone vastly different histories. This suggests that the 


40 


eOpctrollings factor is, not a compositional variation of .the 
erustalisection bur. rather ay Variation “of the physical 
parameters that affect seismic velocities.Whether this is in 
the form of a high temperature gradient regime or in the 
form of a decrease of the effective pressure because of an 
increase of the pore pressure exerted by pore fluids in the 


rock remains problematic. 


4.4 The crust-mantle boundary 

Until rather recently, the Mohorovicic 
discontinuity,the crust mantle boundary, was generally 
considered to be a first order discontinuity where the 
compressional wave velocity increases discontinuously from a 
value of around 7 km/s to a value between 7.8 and 8.6 km/s. 


Although this simple model is still widely used, studies 


with Synthetic seismograms (Fuchs, 1969)combined with 
observations of near vertical reflections from the 
Moho(Clowes,Kanasewich and Cumming, 1969) Suggest a 


Gradational change in velocity rather than a first order 
discontinuity.To explain the observed amplitudes,the change 
in velocity must be accomplished within less than 1km. An 
even more complicated model has been proposed by Clowes and 
Kanasewich (1970) and Meissner(1973) who suggest that the 
crust/mantle transition zone is a many layered complex of 
alternating high and low-velocity laminae having a _ total 
thickness of a few kilometers.Such a sandwich-like structure 


would explain satisfactorily the composite nature of the 
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ret vections trom the moho, that” are’ “often identified in 
field records. It can also explain the anisotropy as well as 


the amplitudes of the Pn waves. 


4.5 Theoretical travel times for a multilayered medium 

A seismic head wave is a diffracted wave,propagating 
along an interface between a lower and a higher constant 
velocity layer.It can often be used to determine the depth 
of the interface as well as the poeta of propagation of 
the elastic waves along the higher velocity layer. The first 
step in solving the problem of a multilayered case is to 
attempt the forward problem.For a completely defined model 
it is a straightforward procedure to compute the theoretical 
travel times of the seismic rays travelling through the 
medium.With the assumption of plane layering and velocity 
increasing with depth one can _ use the travel time 
observations from field records ina direct inversion to 
determine a preliminary model. The travel time calculations 
are approached by applying geometrical optics theory to the 
seismic rays,assuming a multilayered earth in which the 
elastic wave velocity increases monotonically with depth, in 
a stepwise fashion and remaining constant within each layer: 

Case soo <a, (1) 

n is the number of layers including the half space.The free 
Slirface interftacesis, the horizontal” susface upon which the 
source and receivers are situated.The layers may have 


appear ivecioSssrer anuupalp source, (the travel time of the 
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n-layered arbitrarily dipping model is; 


xsinD, n=1L ,-ZU 
Ton = aeons se oe ) on [cosD 


SOS aa OL) eee (25) 
1 G23) Be A Dn 


se Un 


This will give the travel time for a head wave travelling at 
thew toOpTometne = N*th--layers “For the downdip «Source the 
equation is (Adachi,1954): 


xSinU Dede eZ 
qjeors:. Basnl - - 
Dn fe 


LCOS P nCeS Uh. = buon aD 
1 Pal “2 % ‘ eB 


Dn C3) 
Zul and Zdl are the updip and downdip thicknesses of 
the l'th layer respectively.The updip and downdip rays 


,incident on the 1+1 interface, are defined by the angles Dl 


anduu,Ul wathpsrespect..60 .theysvertical, directaon. whrie 


Qn ho ee ela corresponding epeice With pespects to 
the vertical,,.touthe interfaces. 
Dy = Sie + Soa (4) 
By Wa Ron Se pti rol 
6 


where Oe leeisatheadinvangle of «the “laijeth sinterface.The 


Coie ica | retraction sat intertaceun sOocCcurs, when 


at al = Se any es 
“n-1 “n ne 1 
hence 6 Bema nerocner angles Ofsincidence tare scalculated 
n-1 n 


trom -them wise cL Snel.’ s #law of sefraction. The calculation 


starts from ores a and proceeds upwards. 
sind , aoa EEE GD a 3) 
n=? ; Wa-1 
; 7 
sinu _, rae | sin(U_ Pe eae 
a a 
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sind, sin(D,-8.) 
a, a, 
: 8 
sinu, o : sin(U,+8,) 
ae a2 


Equations 2 thru 8 have been programmed for a digital 
computer using FORTRAN language.The computation proceeds 
recursively starting from n=2, until the last layer of the 
model is encountered.The case of a Source burried at some 
depth within the first layer is also taken into account by 
the program. The program will be useful later to test the 
accuracy of the derived models. .Figure 20 shows a ray 
diagram for a three layer model illustrating the geometrical 
Parameters involved in the computations. From eq. (2) and 
(syyptthe slopes of the travel times curves are bin and bh. 


These are the inverse of the apparent velocities. Clearly 


pees oP 1 
b = <a ana b PLS, ye 8! 
D n n 


With field observations the process is rather in 


reverse order.From the apparent velocities one can find 
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The other angles of incidence can be found now using 
equations 7,6 eceim but “in reverse Order “starting from 


1,n, ¥7.n and proceeding downwards.Then from eq. 6,7,4,5 it 
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PovLound) —thatetewscritical angle. (of incidence, at the 


interface n and its dip ® | are given Clore nes )) by: 
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Ray diagram of a three layer model. 
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1 -1(° -1(°n-1 
oa h “we ean e sin(é._5 Vie sin sin(v_ | 
n-2 : n-2 
(1-15) 
rae Bee } sin? [7-3 sin(é jest: ete sin(v ) 
n Oe n-2,n a9 n-2,n 
For n=2 the equations are; (12) 
a a 
64 ora = [sin [+] San =] (lla) 
g ay U2 
a ie a fe 
Qs sin" [54] sin | (12a) 
DZ UZ 


Pe cee ee bye themcirect ) wave-The \velocity (an thes enth 


layer is obtained by Snell's law: 


“n-1 
eo = = CAS) 
n Sen 


The thickness of the n-1 layer is obtained from _ the 


intercept times Idn and Iun and the relations (2) and (3) 


n-2 ees 
: e f 
dee i A (cosD,+cosv,) 


L=1 x 


z. = (14) 
- osD tT COS V 
D,n-1 Cc lb) eT 
n-2 Z 
a Teh) Dis (COSDE COS Y,,)) 
Th ea Tle i a. EN ESF 
= L=1 g 
Z a Nal ty sl a ec ene sel ete eons = 
Viti cosD 7 COSV 
n-l n-l 


Equations 9 thru 15 were also programmed in a_= recursive 
manner to solve the problem of the multilayer dipping 
model.The program will receive as input, information from 
the actual seismic records (apparent velocities and 
intercept times) and will give the depths,dip angles and 


frie velocttpes*eotettic,various layers. it ts interesting to 
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notice that with one source alone,there is insufficient 
information to solve the problem of the n-layer case, there 


beingec\NaW unknowns and only «2(n-1) supnined. by, the data. 


4.6 The t?-x? method 

In. ya wlayer. over..a whali »~space,-the arrival time of 
reflected energy depends not only on the reflector's depth h 
and the velocity v_ above the reLlector but also on ibhe 
offset distance x.The equation,for 5 Single layer over a 


NensGipping “halfispace,relating the,arrival time t,iwwith,x pv 


1S: 
t 2 =s 2 Vee ve (16) 
If the reflecting interface is dipping by an angle 8 an 
extra term shxsint should be added to the right side 
Vv 


of the equation.However for small dip angles (less than 1°) 
and for small offsets this term becomes negligible. Hence, 
if we plot t? versus x? we get a straight line of slope 1/v? 
from which we can determine the velocity v above the 
remlectonshrromethestintercept)ewiths thers tfo tmxdasehwe)iican 
eases the corresponding depth of the interface. For an l 
layered sequence above the reflection it can be shown that 
the velocity v is replaced with a root mean square velocity 
v for the layered sequence. 

Once the RMS velocities v1 and v2 have been determined 
for two successive reflectors, the interval velocity Vi can 
be®found from*Dix (1955) °formula’ 
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wnewe tiyu2-sare® the ’oimnberecepe-' Erimes*26fo@thet esuccessive 
reflection branches. The method assumes straight ray paths 
in the overburden and the existence of a first order 
discontinuity as a reflecting horizon. 

PorPan=accurate’ “calculatYron® of Othe sedepths stand Pethe 
interval velocities one has to confine oneself to the 
subcritical region of the reflection branches. Supercritical 
reflections are affected by refraction and will give a 
larger depth for the reflector as well as larger velocity. 
ALSO, > ar Piangeguobisets¥ ithe eelrech! of Vedipssmaye not be 
negligible.Nevertheless, if the Ful) length of the 
retrograde branch iS used -uSing both subcritical and 
Supercritical reflections then a reduction of 10% for the 
calculated depths is recommended by Giese(1976) and a 
correction of 3% for the average velocities may be applied 
as well.For Subcritical reflections,no correction i15 
needed.As we will see in the interpretation of the data,the 
Te=x2 method proved to be of considerable help in 
determining the depths of reflectors which do not show up as 
refractors and also in revealing information about velocity 


INVEFSLTLONS within -the crust. 


4.7 Synthetic seismograms 

While the general velocity structure of the earth's 
crust can be inferred from studies of the travel times of 
Seismic waves, detailed velocity structure including 
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recognition of low velocity zones, can best be determined by 
a combined interpretation of the travel times and amplitudes 
of refracted and reflected waves. 

Several methods have been used recently to compute 
theoretical seismograms and they have been so far developed 
that synthetic body wave seismograms can be calculated for 
models which are as complicated as those estimated from 
modern travel-time inversions.The main methods are the 
reflectivity method(Fuchs and Mueller,1971),the asymptotic 
ray theory(Hron and Kanasewich,1971) the generalized ray 
method(Wiggins and Helmerger,1974), disc ray theory (Wiggins 
and Madrid, 1974,Wiggins, 1977) and the WKBJ 
seismogram, (Chapman, 1978, Dey-Sarkar and Chapman 1978) which 
visi wsed umy tha'stiwork 

Chapman's method requires a model consisting of 
plane,laterally homogeneous’ and isotropic layers (or 
Spherical shells) overlying a half space. Vertical velocity 
Gradients are approximated by a large number of thin 
homogeneous layers having small velocity contrasts. Layers 
are considered to be perfectly elastic SO that no 
attenuation effects are considered.Shear wave velocities 
were calculated assuming a Poisson's ratio of 0.25 and 
densities are computed fromp=0.252+0.3788V, (Birch, 1954). 
Synthetic vertical component seismograms were computed at 10 
km distance intervals and displayed in record sections using 
asreducinoivelloctty of 62:5 km/sec.Synthetic. and observed 


record sections are plotted on the same scale so that by 
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Superimposing one section on the other, detailed 
characteristics of the observed sections can be compared 
with the theoretical seismograms as a means of 


interpretation. 


4.8 Ray tracing 

The physical basis for ray tracing and travel time 
calculation BOL arbitrary two-dimensional models is 
Huyghen's principle and Snell's law.The following simple 
method of ray tracing was applied in this work: First, a 
convenient system of carteSian coordinates is chosen and the 
equations of the interfaces of the model are defined. 
Then,starting from the source of the seismic energy and a 
given ray parameter,the equation of the ray incident at the 
top of the second layer is derermined.. By successively 
applying Snell's law for this ray, the equations of its 
refracted and reflected segments up to the desired depth are 
defined unital the ray returns to the earth's 
surface.Subsequently, the system of linear equations for the 
ray segments and of the interfaces it has encountered are 
Solved to find the coordinates Of the points of 
intersection. Once these coordinates are known it is easy to 
compute the travel time of the ray by summing up the travel 
time increments amd also to produce a ray diagram uSing a 
computer graphics system. The computation is repeated for a 
large mumber of ray parameters and for various levels of 


maximum penetration. Positive or negative velocity gradients 
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can be satisfactorily approximated by an arbitrary number of 
thin constant velocity layers.Interfaces can be arbitrarily 
dippiug Tand possibility of “wlow, velocity zones, "1S. not 
excluded.The program also computes the times and offsets of 


Ene CHI tcical tellections. 


4.9 Correlation 

Correlation is the first step and perhaps the most 
important one in the interpretation of seismic data. The 
term correlation in seismic interpretation implies the 
identification of the same phase on different seismic traces 
in the time-distance domain. Such an identification of the 
seismic phases is based on the following criteria. 

1) Coherency of the events: arrivals of elastic energy 

from the same seismic phase will have a distinctive 

Similarity in appearance from one recording location to 

another. Clearly, the amplitudes must decrease with 

distance and the shape of the wavelet will change due to 

preferential absorption of higher frequencies at larger 

offsets. 

2) The apparent velocities must show values within a 

physically realizable range. 

3) The travel time branches must be of sufficient length 

to be identified uniquely. 

leo may. be noticed that the criteria for correlation are 
different in reflection and refraction seismology. In 


reflection seismology, correlation of the reflected phases 
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on the time sections is very directly related to the depth 
Om the -crerlector. This. ys not as direct in the 
interpretation of refracted rays. Phase correlation is most 
successful if closely spaced detectors, are used and the 
Spacing between them being usually less than 100 min 
reflection studies. In retraction studies however,the 
spacing of geophones is often a few kilometers. 

TOVDYECkKouUp the @Lirst arrival ) of wave «energy on ‘our 
records we used large scale plots (two inches per second, 
2.5 inches for peak-to-peak large amplitudes). The records 
were not filtered. (Filtering may disturb considerably the 
onset time of the first arrivals). Generally the background 
noise level was low andthe first breaks could easily be 


read with an uncertainty smaller than +0.05 sec. 
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5. RESULTS AND INTERPRETATION 


5.1 Interpretation of line A 

Low noise levels and clear first arrivals events permit 
a quite accurate correlation of the initial waves from both 
the north and the south shot sources of line A.Five linear 
travel time segments were interpreted for each shot.The 
method of least squares was applied to determine their 
Slopes (which are the inverse of the apparent velocities) 
and their intercept and reciprocal times.It is then 
noticed,comparing the two sections, that considerable 
differences in intercept times and apparent velocities 
exist.The intercept times and apparent velocities of the 
northern source are considerably lower than those of the 
southern source.The two sets of seismic observations for 
line A are shown in figure 21 and figure 22 .The travel 
times are plotted with a reduction velocity of 6.5 km/s. 

It should be noted that the average velocity for the 
sedimentary layer using direct waves was observed at only 
one location.Calculation of. v(i) from this single arrival 
time would be highly unreliable because there is a rapid 
increase of seismic velocity through the Mesozoic and 
Paleozoic sedimentary rocks causing significant curvature on 
the travel time branch of the first arriving refracted waves 
to be recorded.Instead,a pseudoaverage v(1) velocity was 
computed using the drill-hole information on depths to the 


Precambrian crystalline basement. AS a consequence of 
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Figure 21.... Record section for line A, shot 1 
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Figure 22.... Record section for shot 8. 
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reciprocity , the extrapolated time of any ray travelling 
from shot point to shot point must be the same when shooting 
in either direction. As can be seen in table 2, good 
agreement between reciprocal times is observed for line A. 
Small corrections to the first arrival times, were applied 
to make the reciprocal times agree and to avoid any 
inconsistencies in the inversion. These time corrections 
were always smaller than the error limits of the first 
arrival times(table 2, columns 8 and 9). The available set 
of first arrivals data of line A is reliable and is suitable 
for the first inversion scheme described in section 4.5. The 
result of this first inversion showing depths,dips and model 
velocities of the layers is shown in table 2. Subsequently, 
the programmed Adachi's equations are used to calculate the 
travel times of the head waves of the derived model and 
compare them with the observed times. It can be seen (table 
29 that agreement is very good and that the only notable 
discrepancy is between the observed and calculated head wave 
times from the basement below the south shot. It is clear 
Ghatesthes Ccalculatedmearrivalistof #itthis *event- assuming 
v(1)=3.30 km/s - are delayed by about 0.1S with respect to 
the observed times. The discrepancy could be removed either 
by reducing the depth of the basement or by in creasing the 
average velocity v(1) below shot 8 .Since a depth of the 
Bbasenent “ato mo akm) 2s (correct “according to drilled well 
information we must accept the second alternative and by 


modelling we find that the average 
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TABLE TZ 


First Inversion Model For Line A 


(V(1)=3.3km/sec) 


The headers of the columns are defined as follows: 


Pe-sinterface. 

S - Shot number. 

N - Number of records used in least squares. 
Rt, - Record t, (sec). 


Mt, — Model t, (sec). 

Depth - Depth of the interface (km). 
V app - Apparent velocity (km/sec). 
RRT - Record Reciprocal Time (sec). 
MRT - Model Reciprocal Time (sec). 


Vtr — True velocity (km/sec). 
DIP - Dip of interface (degrees). 
EN Rt, Mt, alba NL cele, RRT 
Wes 6220/03 70.65 1:3205) 5.822004 50.04 
2 
Sasa Sa20'03; 1:40 302071" 5.93+0'04 SOnit 
iecmee Os t0.04 0.1 O03 4.9+0.3 6.14+0.04 47.91 
3 
S75 1.472004 1:55 3.940.3 6.20+0.04 47.32 
TS CS O108). nS.1 11.520.5 6.41+0.04 46.71 
4 


84 2414006 243 12.5+0.5 6.49+0.04 46.793 


14 2642008 264 18.620.6 6.62+0.04 46.16 
8 4 3.452008 3.46 21:140.6 6.7320.04 46.25 


83 990+014 998 47.8+1.5 840+0.07 44.00 


18 7.4340.10 7.44 36.521.0 7.94+0.04 43.77 
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sedimentary velocity under shot 8 must be 3.5 km/s, in 
contrast with Seo OKM/S = Deloye the North snot: AS willie pe 
seen later, this conclusion is in agreement with the results 
of the interpretation of line C which also suggests a higher 
average sedimentary velocity at the south end of the 
profile. The final -result of our first inversion model , 
based entirely on head waves arrival times, is shown in 
figure 23, illustrating a model of the earth's crust below 
line A, consisting of five layers over a half space. 

Clearly , this model quite satisfactorily reproduces 
the observed first arrival times but itS uniqueness is 
questionable. For instance, low velocity zones, if present 
in the section , cannot be detected by analysis of head 
waves alone since no head wave is possible from the top of 
the low velocity layer. Hence, further modelling , analysis 
of secondary arrivals , ray tracing and synthetic seismogram 
Simulation were carried out in order to investigate the 
uniqueness of the initial model . 

Quite prominent bursts of energy occur at various 
places on both record sections of line A . These secondary 
arrivals were correlated and analysed as reflected energy 
using the T?-X? method as described in section 4.6 to obtain 
depths of the reflectors and interval velocities using Dix's 
formula. It must be mentioned that even after application of 
the 5-17 Hz bandpass filter , the reflections were partly 
Cpseured in ther record sections. A filter of 6-14 Hz made 


the reflections much clearer to observe. The reflections 
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from the basement could not be identified in the records 
GxGeur .Ohlvyeson  theisfirst traces, of -~shot.s]) Gat tiets). 
However, the next two reflection branches (R3 and R4 in 
GUGUTES se 2h n22ee) ware quite, sprominent , ain both asections 


I 


Within the first 70 km. The depths of the reflectors , 
SGCCOnGLnGeetouthewl4—%2 methodseare vat. 5.3 gands 11.6 akmeabelow 
Shoteleanc 5.9mands 12.6 -km sbelows. shot. 8... The interval 
VelOGityeelSwo. Ocukm/s siog.Shotal andi6-07,ekma/s for ushot.é. 
The next step was to compute the synthetic seismograms for 
these first three layers of our model. Agreement between 
observed and synthetic travel times and amplitudes of these 
two reflection branches iS excellent for the following 
depths of the interfaces: 4.9 and 11.6 km below shot 1. and 
Crepelidns hc. /eukMebelowyshoteS pawith ~vic2)= 5.87 km/spnav(3i= 
6.17 km/s. These valueS were incorporated in the _ final 
model. 

Continuing our interpretaion we notice that the fourth 
travel time linear branch in both sections is considerably 
delayed with respect to the previous linear segment and as a 
result of this, a notable break exists in the continuity of 
the first arrivals travel times. On the.other hand,.the weak 
phases observed in both sections at offsets between 110 and 
130km can be interpreted as reflections from an interface at 
16.0 km, the corresponding head wave travel time branch from 
this interface being absent. These features provide a good 
hint on the existence of a low velocity layer in the upper 


crust below line A. Careful examination of the seismograms 
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Shows that,in fact,the weak phases are the continuation of a 
quite prominent reflection branch,R5, clearly observed at 
shorter distances. The reflection branch corresponding to 
the delayed head wave arrivals is also identified and these 
two reflection branches are analyzed using the T?-X? method. 
The results are: depths of reflectors at 16.3 km and 19.9 km 
bedow:, Shot ITand’1728 tkm; 120.41 kmes below! oshot 8 .cetTnterval 
velocities are v(4)=6.33 km/s and v(5)=6.10 km/s which shows 
the presence of a thin low velocity layer situated on the 
topigef the f£itth layer of ourminitialomodel. "Comparison of 
the synthetic and observed reflection branches at the 
boundaries of the low velocity zone indicates some 
disagreement concerning the amplitudes of the near’ vertical 
reflections: Synthetic amplitudes are much greater than the 
observed ones. To remove this discrepancy we decided to make 
the velocity contrasts at the top and the bottom of the low 
velocity zone gradational rather than Airset order 
discontinuities. This has negligible effect on the wide 
angle Be nvene iene amplitudes while reducing considerably the 
amplitudes of the near vertical reflections. The depths of 
the low velocity layer and its interval velocity as 
determined by the synthetic modelling are shown in fig 4 of 
the final model. The layer is slightly dipping to the south 
(dip=0.4°). 

€Continuingeour smodellingtmtowards vgoreatervedepthsy we 
notice quite prominent supercritical reflections, observable 
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reflections from the bottom of the LVZ and earlier than the 
Pm beene te In] Pact Po®8correlation of Cthese events can ibe 
continued without great difficulty toward shorter offsets, 
establishing two new reflection branches R7 and R8 for each 
record (section "of sLine”® A} No head waves arrivals 
corresponding to these reflection branches can be 
identified. The T’?-X? method yields the following depths of 
thee sreflectorss* 26 “km?-+[28¢9@kmebelow *shot = 1c%andm29 15 ukm? 
33.5 km below shot 8. The interval velocity of the layer 
bounded by ghed two interfaces is 6.03km/s indicating a 
second low velocity below line A, this time in the lower 
crust. This layer is embedded within the fifth layer of our 
initial ™*modely which” thas’ *anSsintérvalleveloettymeof.s 6767 
km/s.Synthetic seismogram amplitudes however suggest that 
these velocity contrasts, 6 km/s for LVZ2 and 6.67 km/s 
above and below the layer, are probably incorrect as they 
give rise to very large amplitudes at near vertical 
reflections from the top and the bottom of the low velocity 
zone. On the other hand, our ray tracing calculations 
indicate that a velocity greater than 6.0 km/s would not fit 
the observed times for the reflections at the bottom of the 
low velocity layer. One then would be tempted to increase 
the dip of the base of the LVZ2, which would give good fit 
for one side of the profile but not for the other. Hence, we 
fASTSE OnVthetvelocity ofe6.0km/s for? LVZ2 and we choose to 
fake Mehew@tollowing Ewoe stepsvrinesordery tof reduces, the 


synthetic ‘amplitudes of ‘the near vertical reflections: 
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(PSUDSE I tuseltheafirst onderavellocityadascontinurty wate the 
top of LVZ2 by a gradational velocity change (as in LV21). 
edernecducemenervelocitygcontrastpati the pbottomigot, *LVZ24) by 
considering a velocity of 6.5 km/s below the zone instead of 
6.67 km/s and introduce an average velocity gradient of 0.05 
km/s/km in the lower crust,(0.045 under the north shot and 
0.055 under the south shot). The result of these 
modifications is to obtain acceptable reflection amplitudes 
for LVZ2 and also to reduce considerably the synthetic 
amplitudes of near vertical reflections from the "Moho" 
which otherwise would be quite large in disagreement with 
the ovserved amplitudes. 

Concluding our Peeeone cattion of line A we model the 
reflectionstelrompethesaMohorovicicidiscontinvityl eit ’mustube 
mentioned that the Pn arrivals were of very small amplitude 
but due to the practically zero noise level they were quite 
unambiguously correlated, giving quite accurate intercept 
times and apparent velocities. Hence the results of the 
initial model showing that the Moho is dipping to the south 
by about 2.0 degrees must be reliable. However, the 
calculated depths must be slightly greater than the correct 
ones because of the absence of low velocity layers in the 
initial model. 

Goncelarione Ofemthes (bm. arnavalsersrduiveuhisin both 
record sections of the profile. Unfortunately, because of 
the malfunction of some Sager recording instruments, we 
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Critical region of the Pm branch. However we were able to 
pick the reflections from the Moho mainly guided by rather 
Glear  Ssupercritical= and “subcritical ‘reflections. 29t isvof 
interest at this point, in order to model these reflections, 
BOM cons tden@ethewiresults chmethe ray ‘tracing 20 take inte 
proper account the true attitudes of the interfaces and 
especially the relatively steep dip angle of the moho 
discontinuity. (The synthetic seismogram program considers 
horizontal avertaces Me ein Meact hit Fist Shown’ #tior /rehe 
reflections from the Moho that there is a discrepancy as 
large as 1s between horizontal and true dipping model travel 
times for larger offsets. (At zero offset the two travel 
time branches will coincide). The ray tracing, inluding the 
previously mentioned velocity gradient in the lower crust 
gives depth of the moho at 35.3 km below shot 1 and 46.1 km 
below shot 8. For the purposes of the synthetic amplitude 
calctilation orweae horizontal meMohnoweee ther *interfacey is 
peste vonede@yauesune “depths of the “critically reflecting 
segments of the dipping interface. There 1S a notable 
disagreement between synthetic and observed amplitudes fOr 
the near vertical reflections. In fact, as one can see in 
the seismograms, extremely low amplitudes are observed in 
ene faesesee2s “KnaStromutne shote,foret’ s!e.0se practically 
undistinguishable from noise. These extremely weak near 
vertical reflections can possibly be explained by roughness 
of the moho discontinuity. Also, Synthetic amplitudes for 


*The T?-X? method gives depth of the moho at 38.5 km below 
shot 1 and 48.9 km below shot 8. 


a 
» 

i 

od 
= 


seeds ftehet or 16% 


>) - -_ « » .~ Ge ca 
ne: tae 9 & y Pr 
a « 
r 4 * ~ 
+ es f vt goxe +4 


s j 
fr 4 aw 
« 
a r “u * 
o od 
] ~ 7 
a4 1 Tae 4 
: ; 
° 
, Nw &é 


a hs Decne * “c n£bensie=* 
sa0! - Pa is } 
-_ & a 'y- ; id 
t 
' sone sae uex> my 
7 _ 
a 
se ron. mont sigadey 
a ros 7 
i 7 ~ 7 a it 9 : 7 7 : 
i.pecy ass enol toelIede gs 
<o2 we - : P 
oad 
v 5 “ 7 r 
4 < 4 


64 


the Pm critical region of shot 1 are quite large while no 
Similar high amplitudes are observed in the record section 
of this shot. It iS quite possible however that we would 
have observed them, had the recording at offsets between 120 
and 150 km been successful.figure 24 and figure 25 show the 
final model of line A in cross section and in profile form 
respectively, while figure 26 shows the ray diagram for shot 
1; figure 27 shows the travel time reflection branches as 
calculated by the ray tracing program;figure 28 and figure 


29 show the synthetic seismograms for the final model. 


Two alternative solutions 

For modelling gravity data and for purposes of comparison 
with other continental crustal sections, it is convenient to 
have a simple one or two layer model of the crust. This 
model with a constant velocity or a simple gradient in the 
crust can be used to describe an area which laterally is 
very heterogeneous. | 

ij Starting, from ourmfiamal’ model fon ine A we would 
like to compute an alternative equivalent model (A1) 
Comcustana only of two layers over “the half space: The 
sedimentary cover and an average structure for the crust 
with constant velocity. We invert the reversed Pn arrivals 
assuming a one layer crust to find and average velocity 
(below the sedimentary cover) of 6.4 km/s,or, we can 
ediculate siwerghting”’ factors: fore thesedetailed velocity 


Structure of the final model, in which case we find v crust= 
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Figure 25.... Final model in profile form for line A.The 
numbers near the segments of the velocity7zdepth function are 
the numbers of layers used for the approximation of the 


velocity gradients. 
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Gweocrkiy Sec One Snot) anGuo.5s km, Sec. .Or -Ssiiotec. 
2)Finally we would like to compute a model (A2) of a 
Gonstant velocity Gradientekevathin the crust. The velocity 
function will be V(z)=V,tk(2-z,).. where Zo is the depth of 
the basement and V,=5.87km/sec.Then the one way vertical 


praveiuecimess t) fOr al bayritnavelling through ithisivelocity 


V_+kZ 
Srcuctureswill be ~~ = + 4 Crete where Zmax is the 
Ko Deve eZ 
y We it 1 ' 
deppne ofarhemMonro. yrhen k=" 7 (e0~-1) Ths peguat VOM 1S 
max 


solved by the aid of the computer in an iterative manner to 
find the best value of k that will fit the observed travel 
EimMe Bea 2Maxe 25) /VerUst. erhemnesultimis K=040262 Ssecm™ § Vm 
(=velocity at the vase of the crust)= 6.76 km/sec for shot 1 
Snder=Os02/6 sec Vm = 7.06 km/sec for shot 8.figure 30 


shows the two alternative models in profile form. 


5.2 Interpretation of liné B 

Very low noise level and clear first breaks are 
observed on both record sections of line B which are_ shown 
in figure 31 and figure 32. An exception are the Pn arrivals 
of shot 3 which are severely disturbed by noise.Correlation 
of these events became possible,although still somewhat 
AMBLguUOUS,. alter application “of a Shea fee HZ bandpass 
filcersCorreration Of the first breaks yields five linear 
travel came Segments. 2or shot 4 and ‘six. for —shot. 3. Their 
apparent velocities as well. as intercepts and reciprocal 


times were determined (table 3).It is noticed immediately 
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TABLE 3 


First Inversion Model For Line B. 
(V(1)=3.4km/sec) 


The heacers of the columns are cefined as follows: 


| - Interface. 

S - Shot number. 

N - Number of records used in least squares. 
Rt, — Record t, (sec). 


Mi, - Model t, (sec). 
Depth - Depth of the interface (km). 
V app - Apparent velocity (km/sec). 


RRT -— Record Reciprocal Time (sec). 
MRT - Model Reciprocal Time (sec). 
Vtr -— True velocity (km/sec). 


DIP — Dip of interface (degrees). 
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Figure 30.... Constant velocity and constant velocity 
gradient model for line A 
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Figure -31.... Record section of line B ,snot 3 
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Pigure732.... Record section shot ¢. 
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that unlike line A, agreement of reciprocal times is poor 
except from the first two interfaces and that reciprocal 
times of the SW shot(4) are consistently lower than those of 
the NE shot(3). 

For the purposes of creating a first inversion model 
and because of these discrepancies in the reciprocal times, 
the two sets of data were treated as independent 
sections,assuming that apparent velocities are true 
velocities their differences being entirely due to lateral 
heterogeneity and not to the dips of the interfaces which 
all but the first two will be considered horizontal.The 
first inversion scheme outlined in line A was used to yield 
the depths of the interfaces from the intercept times. The 
result is the initial model shown in fig.1 and table 3. The 
model successfully duplicates the observed head wave times 
but considerable effort was expended in determining the 
extent to which this initial velocity structure satisfied 
other information. 

Reflected phases in line B are quite clear and 
certainly not of the complexity observed in iene A. 
Furthermore, very prominent Crate can reflections are 
observed on the records almost for all the reflection 
branches.These critical reflections provide an additional 
useful constraint for our modelling,since the distances and 
Eimes Of thescritical reflections as” calculated” © for* the 


derived model can be compared with the observed ones. 
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Analysis of the first reflection branch. gives depths of 
the interface at 3.3 km below shot 3 and 4.0 km km _ below 
shot 4,4an wagreement: with «the ‘results: ‘of «the. ‘initial 
model.Interval velocity is 5.70 km/s.Average velocity of the 
sedimentary layer is 3.40 km/s which gives correct depth of 
the basement at both ends of the profile. The second 
reflection branch is clearly observed in both sections 
Micha the dirst 70: km The sicriticale reflection) of» this 
branch Gs tcobsenved -at x=37 km, = k5s:7 for ‘shot 3 echt iis 
not so well defined however for shot 4. Depth of the 
interface is 6.9 km below shot 3 and 9.1 km below shot 4, as 
determined by the ray tracing calculations which also 
Suggest that a slight lateral velocity heterogeneity (6.02 
km/s vs. 6.10 km/s) must be introduced in this layer in 
order to fit satisfactorily the observed times. The 
calculated crestical distance is Xcritical=36 km 
JMierrercal=ts5s for shot ce and Xerirtical=47 km 
PRtcritical—1 296s ‘for tshot a. These \ «figures. sare» ‘hn i egood 
agreement with the observed ones. 

Critical reflections are very prominent for the next 
reflection branch in both record sections.They are observed 
atseterttiv=67 km ef icrdt 2= 190s for shot «3 neand /:Xerre=7in km, 
GVeret 2620S urfomheshot.. 4fneAccardingy tosvthe tray tracing 
travel time calculations the depth of the reflector is at 
£122 km) below «shot. 3 and .13.9 ‘km. -below. shot. 4 The 
calculated values for the ent tcas reflections are: 


Werurt.=626 km>r Bherit=1.90secfor rashot « \3scand’ixerit=20, kn 
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Pe crita2.246— ror “shot §4) “in qood “agreement “with -the 
observed values.It is noticed however that the depth of the 
interface below shot 4, as determined by the reflection 
arrivals-13.9 km- is in notable disagreement with the depth 
as determined by the head wave arrivals -15.8 km-. To 
explain this discrepancy we must first keep in mind that the 
reflecting segments of the interface that give rise to the 
observed retrograde branch (shot 4) are confined within the 
first 50 km ,while the head wave first arrivals information 
comes approximately from the middle of the interface (see 
arrows in figure 33). Hence it is suggested that the 
interface is broken the position of the fault shown in 
figure 34,justified by the fact that no break in the travel 
time branch of the head waves is observed in the _ records. 
The head wave first arrivals (shot 4) are coming from the 
downthrown block of the fault giving rise to a greater 
intercept time which is translated to a depth of 15.8 km. 
The throw of the fault-is about) 1.5 km: Imterval velocity is 
6.23 km/s at the NE side of the profile and 6.30 km/s for 
the SW. Continuing our modelling and after careful 
examination of the seismograms we notice the following. One 
additional reflection branch (R7) arriving before Pm can be 
identified for shot 3. The head wave travel time segment 
corresponding to this branch has an apparent velocity of 
6260.km/s. On the other “hand, two reflection branches RG anc 
R7,arriving before Pm can be identified for shot 4 while 


head wave arrivals corresponding to these reflections are 
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absent in the records. Consequently these two reflection 
branches are analyzed by the T?-X* method which gives depths 
of the interfaces at 26 and 32 km respectively (shot 4). The 
interval velocity was 6.05 km/s while interval velocity 
between R5 and R6 was 6.50 km/s. Hence , R6, R7 are the 
reflections from the top and the bottom of a low velocity 
zone situated in the middle crust below line B. However, 
absence.of R6 for shot 3 shows»|that this) zone is sabsent 
under the NE portion of the profile. This represents a 
Significant lateral heterogeneity in our model. Absence of 
head wave arrivals from the interface at 28.4 km for shot 4 
1s explained by the presence of the LVZ which sufficiently 
delays these events not allowing them to appear as first 
arrivals. Due to their weak amplitudes they could not _ be 
detected as secondary arrivals in our records. 

Concluding our interpretation we would like to map the 
Moho discontinuity of the crust below line B. This was not 
an easy task since reliable information from the Moho was 
limited. Although Pn arrivals for shot 4 were quite clear, 
Pn for shot 3 were very ambiguous. The Pm group was well 
defined for shot 3 but not so well for shot 4 especially for 
the ‘SUbecritacaljwregionwewhere samplztudessewere avery sqlow. 
However, the critical reflection was very clearly observed 
at sXerit=(30 kmpeTycrit=5.40Sa% 

We start from the initial model of a horizontal moho as 
Showneitieche first invession model ;gfige33:—Depth cie-49eckm 


below. shot -3tsand -depth.of 46.5-kmibelow shot 4.,.We,would 
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like now to repeat the initial inversion of the Pn arrivals 
touuineludescthe y effect atot Sethe “hlow bavelocity Uiayer? vewe 
calculate the extra time the rays are spending within the 
low velocity zone to conclude that the depths of the moho 
must be reduced by 0.5 km. Then the new depths will be:48.5 
km and 46 km respectively. 

The next step is to calculate the true dip of the 
interface. We invert the reliable Pn arrivals of shot 4 
assuming VMoho=8.15 km/s. The calculations suggest depth of 
Moho at 40.0 km below shot 3 and 51.7 km below shot 4 (dip 
angle=2.3 °) which gives an excellent agreement between 
calculated and observed Pn arrival times, shot 4. We must 
keep in mind however that these observed Pn arrivals are 
coming from a limited portion of the moho discontinuity. The 
exact position of the refracting segments which give rise to 
the observed first arrivals is shown by the arrows in 
fig.33. Subsequently we would like to model the position of 
the reflecting segments of the interface. Hence we ray trace 
the reflections from the Moho and starting from the SW side 
of the profile we notice a clear disagreement between the 
calculated and observed times. Calculated reflection times 
are significantly greater than observed ones. In particular, 
Pies wcalculared: excritical 16 ati) 142 ckmyand, T critical 2s 
beiUS white, observed kcritical ise U30%km jel scriztical=5.40s 

One might think to remove this discrepancy by increasing 
the velocity v(7). This will eliminate the disagreement 


between calculated and observed critical times. However, it 
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will give an even larger value for Xcritical due to the 
reduced velocity contrast at the base of the crust. 
Therefore, as the only alternative left, we conclude that 
the Moho is faulted below line B, the upthrown block being 
under the SW side of the profile. The throw of the fault is 
about 2.5 km as shown in fig.4 of the final model. Agreement 
between theoretical and observed Pm for the SW side of the 
Mecem 21s NOW Geode Calculated AcritsarssJ75 km, §T’crit: is 
5.55s . Agreement between calculated and observed Pm is also 
satisfactory for the NE side of the model Seiden canere: is 
some indication that the dip angle of the moho at the 
Subcritical region might be less than 2.3°. 

The, tinal, ‘model in’proefide form as shown in figure 35 
;figure 36 is the ray diagram, shot 4,while figure 37 shows 
the travel time calculations;figure 38 and figure 39 show 
the synthetic seismograms for line B, calculated according 
to the final model. Agreement between synthetic and observed 
amplitudes is good. 

Finally, we compute the two alternative models B1 and 
B2 for ame B, following the same steps outlined in the 
interpretation of line A. The results are illustrated in 


figure 40 and table 4. 
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TABLE 4 
Alternative solutions for line B.Depths of Moho and Pn 


velocity same as in model B, 


MODEL SHOT 3 SHOT 4 
Bl Verust=6.44km/sec Verust=6.38km/sec 
B2 K=-0,.03/725ecC >" K=0)0286sec ‘* 


Vmz=7.14km/sec Vn=6.98km/sec 
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Final model B in profile form. 
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Figure 37.... Theoretical travel time Calculations, shot 4. 
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Figure’ 39 «ss. Synthetic “sectron shot 4. 
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Figure 40.... Two alternative solutions for model B. 
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5.3 Interpretation of line C 

The data for line C are presented in figure 41 and 
figure 42.First arrival information was reliable for both 
shots of the profile with the exception of the Pn arrivals 
of the southeast shot (shot 9) which were rather poor. 
Careful correlation of the first breaks yields four linear 
travel time segments for each record section.We start our 
initial inversion assuming an average sedimentary velocity 
Oto aK, Sa OLsoOtnesides Of achesprotile. This gives’ depth 
of the basement at 2.3 km below shot 6 and 2.5 km below shot 
9. However, this latter value is slightly low according to 
basement wells information, which suggests a depth close to 
3 km. This requires an average sedimentary VELOC tye Of, © 3.00 
km/s for the SE side of the profile. Then, calculated head 
wave arrivals for this first layer of our model successfully 
duplicate the corresponding observed times (table 5).The 
next set of observed head wave arrivals observed between 30 
and 90 km is of a split-nature consisting of two parallel to 
each other linear travel time segments, displaced with 
respect .<to .j6éach, sother) -by O,js..<8Thesesdiscontinuities in 
time could be explained neither by a sudden change in 
elevation at the corresponding recording sites nor by near 


surface velocity irregularities since no 
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TABLE #5 


First Inversion Model For Line C. 

(V(1)=3.3km/sec under shot 6) 

(V(1)=3.6km/sec under shot Q) 
The headers of the columns are defined as follows: 


| — Interface. 
S - Shot number. 
N - Number of records used in least squares. 


Ri, —) Recerd it, (see) 

Mi, - Model t, (sec). 

Depth - Depth of the interface (km). 
V app — Apparent velocity (km/sec). 
RRT - Record Reciprocal Time (sec). 
MRT - Model Reciprocal Time (sec). 
Vtr - True velocity (km/sec). 

DIP - Dip of interface (degrees). 
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O30 Siti20:03 si Su tOu5.6 72 C:04ae) 5205) 752130 a 
£2 
62a 420/03 S113 2320. 195:5420 042) S315 45270 
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Figure 42.... Record section shot 6. 
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Similar time displacements are observed for the secondary 
arrivals identified within the same offset range. Hence they 
must be ™hattributed* tosthetrefractorwitserem@and ne fact our 
interpretation is a step-like interface at a depth of 4.5 km 
below shot 9 and 4.0 km below the reverse shot. Interval 
velocity is 5.60 km/s. The observed head wave times from 
this interface, mentioned in table 5, correspond to the near 
the shot arrivals of the split head head wave travel time 
segments. The throw of the steps is less than 1 km. 

The next travel time Segment is very well defined for 
both shots at offsets between 90 and 190 km.ItsS apparent 
vemocuryeiseou4, km/s Eor*shoteS#and’6.62 km/s Tfor +shoe.6 ent 
seems however that these differences of the apparent 
velocities are not associated with the true dip of the 
refractor but rather with differences in true velocities at 
the two sides of the profile. This argument 1S supported by 
the large discrepancies observed in the reciprocal times of 
these arrivals-reciprocal time for the shot 6 travel time 
segment being smaller by 1.1s than that of the reverse shot. 
Hence, for the purposes of a first inversion based on _ head 
wave arrivals we prefer to treat this refractor as 
horizontal and we invert the data assuming true velocity of 
6.62 km/s below the interface, for shot 6, and 6.47 km/s for 
the reverse shot.This yields a depth of 13.4 km which quite 
successfully duplicates the observed intercept and 
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Head wave arrivals from Moho are quite unambiguously 
identified for shot 6, at distances greater than 200 km from 
the shotpoint. Correlation of these events yields an 
apparent velocity of 8.27 km/s and intercept time of 8.9s 
Correlation 1s much more problematic for the reverse. shot 
for which these values can be determined only within a 
rather large margin of uncertainty (see table 5).In spite of 
this however, we still believe that a valid correlation can 
be made*® and we proceed to invert the data assuming an upper 
mantle velocity of 8.10 km/s. This yields depth of the moho 
at 42.4 km below shot 6 and depth of 38.1 km below shot 9, 
GNeseuncertaintyevlin this WatterP@valuembeing t28i5tikm. The 
results of our initial model are fully illustrated in table 
Br 

Quite prominent reflections have been recorded in both 
record sections of line C.The reflection from the basement 
could be clearly identified only in the first record of shot 
9 (offset=6.4 km) at t'=1.4s . This observed time is in good 
agreement with the theoretical value (=1.45s) calculated * 
according? to whe pes ults oflour Whatial, modelwaAnalysis Wot 
the first reflection branch (R3) of the records gives depths 
of the reflector in close agreement with the initially 
derived values. R3 was well defined within the first 40 
km.The next retrograde branch (R4) is of considerable length 
consisting of quite prominent arrivals within the first 130 
ib Rayaeracings caleulatiions’’ show? “that § thet sdepthte Diet the 


>Giving Vapparent=7.94 km/s, t,.=8.145 
*By use of equation (16) Chapter 4. 
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reflector must be somewhat different at the two sides of the 
profile:12.5 km below shot 6 and 14.3 km under the reverse 
SHOUT LACSOTGdnGVtOtOUR Ebnitwallminversion’ “the ‘retiracting 
segments of this interface are positioned horizontally at a 
Geprhpot  i3.4qkmpwhich mayrasuggestwethe spresence snot edsome 
structural relief associated with this interface, the 
Gentraleapartasgohf which ‘smay yaced be more abe less 
horizontal.Agreement between synthetic and observed 
amplitudes , for R4 1S very good for both shots. 
EContinulageRour jinterpretation «Welnotice viormshot 91a 
quite distinctive pattern of reflections arriving 
immediately after the R4 branch.These events are 
particularly strong at distances between 75 and 110 km. 
Furthermore it iS at once noticed that their distinctive 
character sna appearance is strongly similar with the 
reflections associated with the upper low velocity layer of 
line A, (compare, for example, with the arrivals at 
distances between 55 and 90 km of shot 1). Their correlation 
foomntrace | storyeracel<cany imtwkact phbe © icontinuedsiwrrhout 
Qi hic u kty towards shorter offsets, establishing two 
additional retrograde branches R5 and R6.These branches are 
analyzed by the T?-X*? method which gives a low interval 
velocity of 6.11 km/s.Synthetic modelling shows shows that a 
somewhat lower velocity of 6.00 km/s will probably give a 
more adequate agreement between the observed and synthetic 
amp kitudesmatomnathemetop ? ofie thexbVic. akt ois fofminherestrto 


notice that the R5,R6 reflection branches are not observed 


mar i 
oz io aebla ows 
sejevet si7 nebau, 
) Perera 
enizzszie? ed? .noieyeynt 


Pais ms 
° (tt. 


; f ma 
_ 4 —s . : 
»'s * ° , . » 5 
6 ja vllssnosli ted Genétsiesq e168 aa 
j : : / yy 
moz io saises tg ‘et? . t60ebEr 


| 2 t ol a - 
a aale _ 
: Mi To ma? . hs aide © . ie 
5 ’ = be ‘® M 7 - . : 
g ‘ a 3¢ ,* «hd a? ae ae > — mass 
4 \ : > 
<.e 
S08 SSnvinyl anh 2 IPA 


4. onidadidatee (ate ‘o . Tetpode, |) ehzawe2) 


} 
_ Mae el) 
Pave 24 °i? Neswsen rieneetes vat 
: 7 te Je 
selec 4 vs “ inet 


98 


for the NW shot which suggests that the presence of the low 
velocity layer is justified only for the SE side of the 
prifile.Probably it is the continuation of the upper low 
Velocity layer-identified “under line A. 

Reflection arrivals from the Moho are designated by 
quite * strong@tamplitudes ini? the critical’ regrvonefor (both 
record sections of the line.However, before modelling the Pm 
arrivals we would like to focus our attention into a well 
defined set of events appearing in the records of shot 6 at 
distances around #50 kiteface teach SshAecAlthough theis 
amplitudes are low, their coherency 1S adequately good _ so 
that we feel that they should not be neglected in our 
modelling. Their eer eieccn from trace to trace establishes 
the R7 reflection branch shown in fig.42.No corresponding 
head waves are identified in the records.For the purposes of 
our ray tracing modelling we consider the interface as 
MROeuzontal Since Nowmsurticient scontrol. of) the, (dip. lis 
provided by the data.The calculations suggest a depth of 31 
km and an interval velocity between the R4 and R7 reflectors 
close to 6.58 km/s. The interface is absent under the SE 
Side of the profile since no corresponding reflections or 
head waves are observed for shot 9. Synthetic amplitudes for 
R7 are in agreement with the observed ones. 

Concluding our interpretation we model the Pm arrivals 
observed in the records.According to the ray bracing 
calculations the depth of the Moho is at 39.5 km below shot 


9 and 44 km below shot 6.Seismic velocity at the base of the 
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Grust 15 6.85 km/sec below shot 6 and 6.65 km/sec below the 
Be verce Shot-Agqreement between synthetic Pm amplitudes “and 
Ceeeavea. SONGS 81S) SatIStaCtOry a es OMe can see in the 
Synthetic Sections firgurené3 jand @voqures44 the final model C 
1s shown in figure 45. 

One —cannot overlook the-spresence s-of-gquite strong 
@erivalStlappeawing, in the records after the Pm curve.Their 
COurelation Voevery § good ‘establishing several ‘retrograde 
curves after the Pm arrivals.The apparent velocity of these 
abriVvalS iSewwhLoOnywand.increases: with déepth...All these 
branches seem to merge near the critical region giving rise 
to particularly strong amplitudes and consequently are 
VOoatnvyearcorerurt ner wa tentvron aACccoraMmatomtnerresutts of ‘our 
modelling they are interpreted as due to a laminated 
SUD-Mond;, Structure €onsisting of Several high }and low 
velocity layers sandwiched between each other, as shown in 
the final.—profate—_.model) figure._£.46..Synthetic...amplitudes 
computed fot this sub-moho stucture are shown for shot 9 in 
igure: 4), while unew COlrresponaing ray Calculations are 
included in figure 48 

The two alternative models C1 and C2 are illustrated in 
table 6 and in figure 49. TABLE 6 

Alternative solutions for line C.Depths of Moho and Pn 
velocity same as in model C. 


MODEL SHOT 6 SHOT 9 
C4 Verust=6.52km/sec Verust=6.41km/sec 
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Bocure 25... bitedemode! for lines C.Solic lines indicate 
where head waves provide direct evidence for an interface. 
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Figure 47.... Synthetic seismograms,shot OM nodelvc 
including a laminated sub-moho structure.Notice the higher 
empiatides of thevem branch in ‘the critical region. 
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6. CONCLUSIONS AND DISCUSSION 

The results of the interpretation of lines A,B,C of the 
1981 CO-CRUST experiment reveal a very complex structure of 
the crust below South-Central Saskatchewan.Low velocity 
layers have been detected below all three profiles.This is a 
particularly interesting feature of the derived models since 
no such layers have been identified by previous seismic 
SurveyS in Alberta,Manitoba and eastern and western 
Saskatchwean.It is of interest to speculate on the nature of 
these low velocity zones and we will return to this matter 
at the end of this section.The crust thickens southwards and 
thins eastwards.Crustal thickness below the study area 
varies from 37 km to 47 km.Depths of Moho below the triangle 
are shown infigure 50.No evidence of anisotropy Or 
Substantial variation of the Pn velocity can be postulated 
Since upper mantle velocity was found to be practically the 
Same below the three profiles where reversed travel times 
were available.The Pn velocity was close to an average value 
of 8.13 km/s.This value is in agreement with the results of 
other selsmic Studies in Alberta and 
Saskatchewan. (Substantially lower upper mantle velocities 
have been reported for Manitoba). This similarity of density 
values of upper mantle material suggests that the base of 
the lithosphere below Western Canada is subjected to a 
uniform thermal regime.There is no evidence of any localized 
partial melting and thermal alteration of the lower crust.It 


must be mentioned that such evidence has been reported 
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previously by seismic studies in U.S.(e.g. in the Basin and 
Range area,Gish et al.,1981) and was explained as a result 
of magma intrusion into the base of the lithosphere caused 
by upward convection within the asthenosphere(Withjack, 
1979) .However,a slight lateral velocity variation is present 
an'@"most= of the layers” of Your +modéls*B and Ci showing™a 
consistent westward increase of eareeata te within the 
crust.This may be due toa slightly elevated temperature 
regime within the crust of eastern Saskatchewan. 

There is clear evidence that the base of the crust is 
faulted under line B.The throw of the fault is at least 3 
km.A Similar structural variation at the intersection of 
lines A and C shows that the fault is probably continued 
through southern Saskatchewan,passing through the earthquake 
prone area of Bengough.The great depth of the Moho near 
Swift current-below shot 4 of line B can be attributed to 
the complex N.E-S.W. trends in Saskatchwean intersecting the 
buried east-west Precambrian rift in Alberta discovered by 
Kanasewich et al.(1969). Correlation between the various 
layers of models A,B and C is reasonable both in depth and 
structure.Considerable effort was spent in attempting to 
derive®©oether? models“? that’ “would fit “the *detaats* of the 
observed seismograms. All these efforts were unsuccessful.A 
total of about 200 models were tested by synthetic 
seismograms and ray tracing before concluding the final 
models appearing in figures 24,34,45. It 1s noted that a 7 


km/s’ layer is absént in-all three models.It 1s also absent 
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Py thesinterpretation of the CO-CRUST C-1979 ‘osrofile which 
is in very good agreement with our model A. 

As mentioned earlier ,line A was located close to the 
NACP conductivity anomaly as postulated by Alabi et 
ate, )2)o,see figure 51. ‘According to their (results “the 
anomaly is a major feature in the lithosphere of Precambrian 
north America, extending from north-central Saskatchewan to 
the Black Hills of Wyoming and thence to the northern 
Southern Rockies.The southern part of the anomaly which is a 
Marrow belt(around 50 km wide) of very high electrical 
conductivity was modelled by Camfield and Gough(1977) to be 
in the upper mantle while its maximum depth in Saskatchewan 
1S not more than 35 km. No anomalously high heat flow has 
been reported for the area of the NACP anomaly in 
Saskatchewan.Therefore the observed high conductivity must 
be attributed to different rock type rather than different 
temperature environment. Whether this is in the form of 
gGraphitic sheets, or brines rich in heavy metals transferred 
iMcOe Tracture fs Zones, FOr uihigh, “concentration of hydrated 
minerals in the host rock is a problem which cannot be 
resolved by the present amount of information. 

AcCcCOrdinGa-comsour—-worktwor=low==velocity < zones: were 
identified in the Earth's crust below line A, at depths of 
15 and 25 km respectively.It is possible that both of these 
low velocity layers contribute to the observed conductivity 
anomaly.There is some indication (see final model C) that 


the conductor extends itself under the SE side of line C. 
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Ponce toe.) bOcatwon. of the conductive body in noreh 
american Central plains,after Alabi,Camfield and Gough 1975 


ja fe 


Camfield and Gough (1977) postulated that the linear crustal 
Structure may be a major continental fracture zone mapped 
Over 2 totale length Of 1000 Km lt ise ing tact quitespossibie 
that the fault we have detected in the crust under the 
triangle ABC(its strike being in alignment with the NACP 


body) comes in support of this argument. 
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APPENDIX 


LISTING OFA COMPUTER PROGRAMS 
Gespiking program 
POrwalcmsowution for ea n-layered system. 
invert Semsotution for a mo=leayered system. 


ray trecang.s program 


ile. 
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THE DESPIKING PROGRAM 
ENTEGERS <2 -1T (7 6B LO oe el 1 Lo 
Crimerhe Ol Dek ee OMe hen 
POL )} 

LBZ = 2") 
PGRRe1S1GN4 15401} 
IGR=ISIGN(1,1D2) 
Gas 


J-W9 


PEL DMA EER 
EPO EA EE 
Pr eGd sinew 0) GOrTOs60 
IT(U)=IGS*IRd 

GO 16:8 
PBEZ= ED id Si d- 1) 
LDEZSIABSHIDEZ ) 

PE ( DeZ= i pet Rd} GO. i 0e22 
HewiGie b1.0)GO 10 24 
Pid) =265*511 

GOV iG" :2 

CONTINUE 

IT(J)=I1D(J) 

IGRR=IGR 

IGR=IGS 

ITS=IT(u) 
PGS-1TS51GN(1,1TS) 
1G3=1G2 

IG2=I1G1 

1G1=1G0 

1GO=1 

Gis tOG3s 

CONTINUE 

LIS 1GS*511 

IGRR=IGR 

IGR=1GS 

DRSeit (uy 

PGSs tS EGNETS TTS) 
IG2=I1G1 

1Gi=1G0 

1GO=(-1) 

GON TO B23 

CONTINUE 
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o wooo 
oa a = 


> are ¢ 


903 
504 


505 
5900 


600 


550 
50 1 
502 


IF(1G2.GT.0)GO TO 23 
IT(U)=IGRR*(-511) 

GO TO 8 

CONTINUE 

Ci Guen Diu: 

GO TO 9 

CONTINUE 

Di Gus GR - 5:14} 

GO TO 8 

CONTINUE | 
Teeter On Dia). 11 (es GS= GR. TGRR,'G0., 1G1, 162,763 
FORMAT(2110,716) 

CONTINUE 

CONT INUE 

O@550 00 145766... 


SOs w ateaee GO pO 500 


) 
pay Gal eNe. (5511) ) GO FG 500 
Fegtl ol teen Ge ae EO: S00 
WRITE. 503) 1 Del DCL em eet?) 
WRITE(6,504)1T(1), T(1+1), 17 (1+2) 


FORMAT(’ ID ahaa i TRAP’ , 318) 
FORMAT(’ IT VALUES IN TRAP’ ,318/) 
Tal, alee ta) =o Dal Ser 1.) 
WeiheoeSUoN leach). li (ie tll pis?) 
FORMAT(‘’ MODIFIED IT VALUES ‘ ,318/) 
CONTINUE 
ae 600 K=4,766 
F(IT(K).LT.495)GO TO 600 
ese res -511))G0 To7600 
Feliu ose e501 GO 10 600 
oda Weert 
CONTINUE 
DO sO eal GOL S 
IPP=I1+7 
WRiEd.6,.901)(1D(K) Kel, TPP) 
WRITE(6,502)(IT(K) ,K=1,1PP) 
CONTINUE 
FORMAT (818) 
FORMAT(818/) 
STOP 
END 
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THE FOREWARD SOLUTION OF A MULTILAYERED SYSteEmM. 
THE CALCULATION CAN BE PEFORMED FOR AN ARBITRARILY LARGE NUMBER 
OF LAYERS. THE LAYERS CAN BE ARBITRARILY DIPPING.NO LOW VELOCITY 
LAYERS ARE ALLOWED IN THE MODEL.THE SOURCE CAN BE AT ANY DEPTH 
JITHIN THE FIRST LAYER.THE RESULTS OF THE TRAVEL-TIME 
CALCULATIONS ARE PLOTTED ON A T-X AND A T’-X GRAPH. 
DIMENSION PA Gey), ZU 0} <2), DIP C10). D010),U116)} 
DIMENSION XSA(32),TSA(32),RSA(32) ,ASX(32) 
INTEGER VAR‘ 
READ(5,31)A 
FORMAT (10F5.2) 
READ(5,32)2ZU 
FORMAT(10F6.3) 
READ(5,33)ZD 
FORMAT (10F6.3) 
READ(5,34)DIP 
FORMAT(10F5.2) 
READ(5,110)PL,PH,PARAM1, PARAM2 ,MODEL 
QO FORMAT(4F5.1,14) 
READ(5,120)VAR1 
0 FORMAT(14) 
READ(5,124)VAR3,OFFSET,RVEL 
4 FORMAT(3F7.2) 
MODEL 1S THE NUMBER OF LAYERS OVER THE HALF SPACE 


DUMENSLONARMATO)8 DIUMOUR RNI 10}, 014490), TUL10),1TO(10) RIU(10) ,F 
R) 


CARIB ALOIS 
CALL (PLOmi 
CALL <P LOG 
CALLS PLONE 
CACHE PLOIE 
CALL PLOTA 
CALL pel 
( 
( 
( 
( 


~— (ee) 


— 


GALE ePEGn. 
CALL PULOR 
CALL PLOT 
CALE TREO 
CALL PLOT 
NSY=0 
TPL=1 Petr) 
IPH=IFIX(PH) 
VAR2=290.-OFFSET 

ALFA IS THE DEPTH OF THE SOURCE 
ALEA=0 20 
ZU( t= 2001 Je AURA 
ZD(1 = Z DUA BA 

CONVERT THE ANGLES TO RADIANS 
DO 1 NA=8)10 
DIP(NA)=DIP(NA)*0.01745 

1 CONTINUE 

CALCULATE ANGLES D(10),U(10) 
DO 42 LEK=1,MODEL 
D(LEK)=ARSIN(A(LEK)/A(LEK+1) )+DIP(LEK+1) 
UE =ERoENBAC UE) / NI LEKI) J-DIPCLEK+1) 
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42 CONTINUE 
CeDIREC 1 aWwAVE 


WRITE(6,89) 
83 FORMAT(//,’DIRECT WAVE’ ) 
NCPL=0 


DOBRA =10 -VAR LT, 10 

NCPL=NCPL+ 1 

ket 10: 

ee Wabi 1290.) X=VARS 
=(SQRT((ALFA*ALFA)+(X*X)))/A(1) 
Ae DIDI X/RVEL) 

WRIRE GS 3.00.) 12, RED LO 


XSA(NCPL) =x 
ASX(NCPL)=OFFSET-X 
USARNGCR LA) =F 


RSACNGPE) =REDIDI 
SUC. BORMAL GE CU aar 10). 2) FO anes 4) 
5 CONTINUE 


RMSA=0.0 

TMSA=0.0 

DO Shit. Wee Coke 

BAe Se KyeoG1 . 1MSA)TMSA= TSA CK} 
Et RAG ar. RMSA)RMSA=RSA(K) 


oe 1 GON TA NUE 


RrUMSAL Gi.60)) GO Ye Sy 
T2P=(((IFIX(TMSA) )/2+1)*2) 
GOP TO a2 


Sep eter= (ULF IACIMSA) )/10)41)*10 
go 2 CONTINUE 
IF(RMSA.GT.10)GO TO 54 
R2P=(((1FIX(RMSA))/2)+1)*2 
GOs T0558 
54 R2P=(((IFIX(RMSA))/10)+1)*10 
58 CONTINUE 
K2P=OFP SEW /PL 
XSA(NCPL+1)=0.0 


XSA(NCPL+2)= Aer 
ASX(NCPL+1)= 
ASK(NGPL=2) =x3F 
TSA(NCPL+1)=0 
TSA(NGPL #2) 217 /Ph 
RSA(NCPL+1)=0 


RSA(NCPL+2)= bape 
231 FORMAT(8F10.4) 
CALL LINEP(0.042) 
CAE LANE (XSA ShSde NCP MEWS NSY) 
CALEP TENET ASXS IGA. NCPL.,1,1,NSY) 
DO 4 KA=1,I1PH 
RL=KA*TSA(NCPL+2) 
RKA=KA 
CALL NUMBER(-0.35,RKA,0.086,RL, ae 0,-1) 
GALL SYMBOL(0..0,, RKA, 010425 of ue nt) 
CARLPSYMBOL(PL., RKA, 0. G42 ce te) 
4 CONTINUE 
DOv.322 Waa4 TPL 
RL=(KA*XSA(NCPL+2))/1000. 
RKL=KA 
CALL NUMBER(RKA-0.12,-0.20,0.087,RL,0.0,1) 
CAMPASYNEGICRKAA0T0. 0.042,3,0:0.1 
322 CONTINUE 
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100 
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Chia PleoiguP lot: 5).,.0..0e-- 33) 122 
DO 15 KA=1,1PH 
RL=KA*RSA(NCPL+2) 


RKA=KA 

CAG le sate Oral Oe RK A O04 27 37 me OS, 
Gaile oni CPL RKA AG 042, 2 Digi) 
CAamNOMBe RG -0..35:, RKA Oe 086, eee Oye) 
CONTINUE 


De gat Sen maid PL 

Res KOA UNCPL +2) )/ 1000" 

RKA=KA 

Chane O. 2m On206O.0G7, RL, 02051) 
CALS VNCOERKAY O10) 107, 0425535.07.0 , 1) 
CONTINUE 

Cie NEoweR SA, NCPLest, 1s, 
CANE enon ROA, NCPL. diet, 
CARP eGiiiee Rb t.5 | 00), 734 
NSY=NSY+1 

eee Ar ANeaNerOn.0) GOP 1G) 80 
VAR DABEE SPARAM 1S NOT EQUAL 710 0. OVEN HE SOURCE 


Pome Ae eserieem. SECOND 7 INTERPACE 
URDEP. SOURCE 
Core Cemererr one rO. THE sShDEGOr JHE PROFILE IN WHICH ZU 


DEFINED.(IF ZU>ZD DIP<0.0). 

Do SAR Sia) /At 2) } 

Weine (io. 701) 

FORMAT(//,’UPDIP SOURCE HEAD WAVE FROM INTERFACE 2’ ) 
NCPL=0 

Dio we ho Rd 2 0 

NCPL=NCPL+1 
X=1X-10-(ALFA*TA 
IF(X.EQ.VAR1-10- 
-(ALFA*TAN(D(1) 
FUSS ean eleeS 


) ) 
TAN(D(1))))X=VAR1- 10 
Ui ALI) Ve UcOS PDL T))) 


TUZ=TU2ZS (AGRA, A 
REDM= 1. U2—UXRVE 
TT 1=X+(ALFA*TAN 
WRITE(6,100)TT1 
FORMAT (’ 0! , Ee 
XSA(NCPL) = 1T1 
TSA(NCPL) =TU2 
RSA(NCPL)=REDT 
CONTINUE 
Come a one OAMNCPL 1 
CAEL PLOT(PL+1. oF 0. 02-3) 
al 
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ET 
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ee 


CALERA NEL XSAc RSA, NCP 
Ce eee (PL+1, 5), 0.0 
NSY=NSY+1 


C DOWN DIP SOURCE 
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NCPL=0 

WRITE(6,702) 

FORMAT(//,’DOWNDIP SOURCE HEAD WAVE FROM INTERFACE 2’ ) 
Dog hae Rae 10 

NCPL=NCPL+1 

X=1X-10.-(ALFA*TAN(U(1)}) 

Looe OVER 10 1h) PATON OU 4) RX EVAR 1210 
-(ALFA*TAN(U(1)))-VAR2 

Too ee oN Ut 2001) AL te (COS (Ot t+ 
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EU eye ele) CCOS CUI: Tok) (ea 
Poe 2 SCAR RA GAR GOS (UTI 

TOZ=TDIOS RARE ALAN) STAN (UIT) eSINCU(1)) 

REDTD=TD2-(X/RVEL)-( (ALFA/RVEL)*TAN(U(1))) 
TT1=X+(ALFA*TAN(U(1))) 

WRITE (CO 2200) TT14.. 102) REDID 


ASACNCP I] = 111 

FSAI NCPL )=TD2 

RSA(NCPL ) =REDTD 

FORMAL (Of RwO.2,7F8.4,F8%4) 
CONTINUE 

CAME INECASKM] SAENCPLES a aNSY) 
CAMS SPOT Pie i500, >3) 

GAS LINE CASA RSA, NCPLA Tat NST) 
GAGE PEOTA =. Phen. 5) 8070)-3) 


C UIHRERBEAY ERS 
CoUPDTeSSBURCE 


* 80 


90 


400 
645 


Ue Cee MRA cock RCD 


W238 =NRSIN(CA(2) /A13} ) 

Bi= (ACT) /A(2)) ¥SIN( D1 2))=DIP( 28 } 
Hise ARSIN(B1) 

Cra CGwy hie eS INU U2)SD1P (2+) 
BAS steer auch 


REAL NEWD1 
REAL NEWU1 

NEWD1=D13+DIP(2) 

NEWU1=U13-DIP(2) 

WRITE(6,930) 

FORMAT(//,’UPDIP SOURCE HEAD WAVE FROM INTERFACE 3° ) 
NCPL=0 

POLG4IS01X=102VAR1,10 

NCPL=NCPL+1 

X=1X-10.- (ALFA*TAN(NEWD1) ) 
IF(X.EQ.VAR1-10-(ALFA*TAN(NEWD1)))X=VAR1-10 


mca eee: NEMA! 


F-(ALFA*TAN(NEWD1))-VAR2 


TU3=(X/A(1))*SIN(NEWD1)+(ZU(1 


TU3=TU3+(ALFA/A(1))*COS(NEWD1 
NEWD1) ) 


TU3=TU3+( (ALFA/A(1))*TAN(NEWD 1 
EWD1)) 


)*S] 
RTU3=TU3- (X/RVEL) -((ALFA/RVEL) *TAN 
TT1=X+(ALFA*TAN(NEWD1) ) 
WRITENG, 400) 711.103, RTU3 
KSATNGRES Eta 
TSA(NCPL)=TU3 
RSA(NCPL ) =RTU3 
PORMATING 07 SFO 2ihee 4 FS. 4) 

CONTINUE 
OV aUNE KSA TSANEPL 1. 1,NSY) 
CAC PPLOF RUSS t OHO +3) 
CALE VINE XSA,RSA,NCPL,1,1,NSY) 
ChUEV POTS (PU+1 5)0.0,-3) 
NSY=NSY+1 
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C DOWNDIP SOURCE 
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WRITE(6,91) 

FORMAT(//,’DOWNDIP SOURCE HEAD WAVE FROM INTERFACE 3’ ) 
NCPL=0 

DO’ 7 1X=10;VAR1, 10 

NCPL=NCPL+1 

X=1X-10.-(ALFA*TAN(NEWU1) ) 
IF(X.EQ.VAR1-10-(ALFA*TAN(NEWU1)))X=VAR1-10 
G- (ALFA*TAN(NEWU1))-VAR2 
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TORIES Xa AK 
N+(ZD(2)/A 
1P3= 13+ { 
IDSs TDs+( 
RI Ds=TD3- 


1) 
2 
Ve 
(A 
( 


)=SIN(NEWUt ) + Daag /'A te ee ee eC ENUM) 
)}) 


) ) Urcaeniear eae U(2 


Pe seca nett 
LFA/A(1))*TAN(NEWU1)*SIN(NEWUT1) ) 
X/RVEL)-((ALFA/RVEL)*TAN(NEWU1) ) 


TT1=X+(ALFA*TAN(NEWU1) ) 


WRITE(6,500 


THs ADSBRT DS 


ASAM NG PIE T7 G 
TSA(NCPL)=TD3 
RSA(NCPL)=RTD3 


FORMAT (’ 0’ 
CONTINUE 

CALL LINE 
CABba PLOT 


/F10.2,F8.4,F8.4) 


ASK Bee NICE TINS 1) 
PL+1. oe. 0, -3 


~— 


( 
( 

GAERPLINEIASX. a NCPL, ASG, 
(rt 


COLE IREOT 
NSY=NSY+1 


(Pile i. 5), OF 0, 73)) 


PARAM2 ISNOT EQUAL TO 0.0 WHEN THE MODEL CONSISTS ONLY OF 
POR ER EEAMERS 
TEUPARAM2Z¢NE*0.0)GO TO. 22 

GENERALIZATION FOR ANY NUMBER OF LAYERS 

DO 20 J=3,MODEL 

IN ITS GENERAL FORM,THE PREVIOUS STATEMENT IS 
DO 20 J=3,N WHERE N=(TOTAL NUMB. OF LAYERS-1) 


AQ 


60 
cM 


IVAL=uU+1 
RM(1)=0. 
RM(2) = anotN 
INDEX 1= 
ma eee 
DI(J-INDEX1 
INDEX=1 
DO 40 K=3,u 
RM(K)=ARSIN 
M-DIP(JU-INDE 
INDEX1=INDE 
IF({JU-INDEX1 
DI(U-INDEX1 
INDEX=INDEX 
CONTINUE 
RN(1)=0. 
RN(2)=ARSIN 
IND1= 
elas) = 0. 
UI(J-IND1) 
IND=1 
DOVSO KES. J 
RN(K) = ARSIN 
N+DIP(J-IND) 
Rhee 
TRA USIND PL 
UI{J-IND1)= 
IND=IND+1 
CONTINUE 
WRITE(6,92) 
FORMAT(//,’ 
NCPL=0 


(A(J-1)/A(U)*SIN(D(U)-DIP(J))) 


)=RM(2)+DIP(u) 


ieee ees ew ies TNOEN 
X 


X1+1 
PET oNGe eo) 40 
)=RM(K)+DIP(J-INDEX) 
+4 


(A(U-1)/A(U)*SIN(U(JU)+DIP(u))) 


=RN(2)-DIP(J) 


eam Gh Cru ean ECEGSAINDY 


a ae ineas 16) 
RN(K)-DIP(JU-IND) 


IVAL 
DED SOURCES HEAD WAVE FROM INTERFACE’ 


DOP Sea Kesh ARA BO 


NCPL=NCPL+1 


X=1X-10.-(ALFA*TAN(RM(JU)+DIP(2 


LEEKS BORN OR 


))) 
inet AL EAST ANURM( J )4+DI1P(2))))X%=VAR1=-10 


,14) 


124 
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_ Wu . 4 
{ i] *) 2OAAT AL NORA 7 ss rh or r 


TAN(RM(J)+DIP(2))] 


E-(ALFA* “VARZ 

TUG =O: 

TU(2)=0. 

TU(3)=0. 

K=u+ 1 

PU CAA SIN CRN REP A) 

Th) =o hee EA ACOS CRM Cds aber ee) AT) 
TOO rele CZ) 7 A (1) Je CCOS CRM} +DIP (2) )+COS(RN( J) 
P= DUPA2Ie))) 

TUCK STUCK) +((ALFA/A(1))*TAN(RM(Jd)+D1IP(2))*SIN(RM(U) 
ESO 2 ai) 

TUG = TUCK AU Ata ee (COSI Dd) )+COSTUCd) )) 
PNREGER SR 

R=J-1 

DOBS 2 Sle 2e0k 

Te nei) eaeuime ACT) eVC OSLO I AI) )+COS(UT{1))) 


a2 CONTINUE 
RUCK) SEU eX RIVE ea A 
TT1=X+(ALFA*TAN(RM(JU)+DIP 
XSA(NCPL ) = ie 


TSAUNCEE ) =a.UiKs) 
RSA(NCPL)= ealiie 
WRITE(6,71)171,7U(K) ,RTU( 
71 FORMAT(‘ 0’ ,F10.2,F8.4,F8 
53 CONTINUE 
Ghee X SAM ESAANCeL 1: 
CALL PLOT(PL+1.5,0.0,-3) 
CAME VEINE CX SAS RSA,NCPL, 17, 
CALL PLOT (—tPL+4.5)>..0.0, - 
NSY=NSY+1 


C DOWNDIP SOURCE 
WRITE(6,93)I1VAL 
23 FORMAT(//,‘ DOWNDIP SOURCE 

NCPL=0 : 
DOsos al Xs ORVAR MS LO 
NCPL=NCPL+1 

X=1X-10.- (ALFA*TAN(RN(J) 
DE Xe EO VARS TO =ALE As TA 
Ke tie hl astra 
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1 
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Aes 
eZ Asin) 3) 


ual GO) Ty 00) ESS ene 00 


he 
SWABS ACA ey 
K)+(ALFA*COS 
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62 CONT 
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1 et ROK) sR 
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ott eee 
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ae ee eager) 
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1,NSY) 
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THEAD WAVE FROM INTERFACE: 


BUR ies) 
N(RN(Ju)-DIP(2))))X=VAR1 

2 lea Aree 

-DIP(2)) 
*(COS(RM(JUJ)+DIP(2))+COS(RN(u) 


TAN(RN(uU)-DIP(2))*SIN(RN(J) 
D 


LCOS Poms COS( O02) 


ieee 


S 
D Veet 

LFA/RVEL)*TAN(RN(JU)-DIP(2))) 

ee 
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26 


RSA(NCPL)=RTD(K) 

PO RNieln Os ara. 2. Old, Pe. 4) 
CONTINUE 

CALL LINE(ASX, Wee NG Pierre NN Se) 
eb igoroy oka AO 0, az 

Citta (ALNE CAS Me peu NCPL, shea RNR) Ga, 
Sie ecient, 5), 0. 0, -3) 
NSY=NSY#1 

CONTINUE 

CONTINUE 

GAEL: «PLGuel 03.c0-.0;40,999)) 

STOP 

ENN 
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C THE INVERSE SOLUTION OF A MULTILAYER SYSTEM. 
Ce THE PROGRAMS CAUCULATES DIPS, VELOCITIES AND UPDIP AND 
C .DOWNDIP DEPTHS FROM INTERCEPT TIMES AND APPARENT VELOCITIES 
Gy (NTOTAL 1S THE NUMBER OF LAYERS INCLUDING THE HALF SPACE 
G 
DONE NST ONAL Ole Vid tO RID ULOy eR eet Ole Alen DIP (90) 
DIMENSION D(1 ead [1Od BDMaaed AUN, ZD110) .Z0010) 
DIMENSION DL(10),UL(10), eek 
READ(5,1)DVEL, NTOTAL 
1 FORMAT(F6. aa ey 
Wie = DVE ls 
READ CS. QV(VDCK Ke 1, 7) 
2 FORMAT(7F6.2) 
READ(5,2)(VU(KK) ,KK=1,7) 
PERDS 2 Uk EO ene es lone} 
PEND E522 hur luigels) .LL= te fl 
D12=0.5*(ARSIN(A(1)/VD(2))+ARSIN(A(1)/VU(2) ) )} 
DIP(2)=0.5*(ARSIN(A(1)/VD(2))-ARSIN(A(1)/VU(2))) 
ANGLE (2)=57.2958*DIP(2) 
WRITE(6,3) ANGLE (2) 
OME OR VAT Cw bodes" bo wor sD GGRGES thi 
D(1)=D12+DIP(2) 
UEP = Dito br eaters) 
Nee A) SUN Dh) 
WRITE(6,4)A(2) 
AO BEGr EO! eee y= er oO. 2.” KM SEC 5/y) 
7D(1)=(A(1)*RID(2)}/(COS(D Per +cost(UCt)y) 
OMe eet ae RLU iy CeOS (D114) )+COS(U11))) 
MPEG 6, 51.2 ae 
ROPER RWC 10, oe eae iar ae. 3, KM /) 
WRITE(6,6)ZU(1) 
6 CRORMAN (07%. ZU Gh) = . FS. 88 KN Ys) 
C THREE LAYERS 
D13=(-1)*DIP(2)+ARSIN(A(1)/VD(3)) 
U13=DIP(2)+ARSIN(A(1)/VU(3)) 
DD(1)=D13+DIP(2) 
UU(1)=U13-D1P(2) 
D23=0.5*(ARSIN(A(2)*SIN(D13)/A(1))+ 
FARSIN(A(2)*SIN(U13)/A(1))) 
DIP(3)=DIP(2)+0.5*(ARSIN(A(2)*SIN(D13)/A(1)) 
E-ARSIN(A(2)*SIN(U13)/A(1))) 
ANGLE (3)=57.2958*DIP(3) 
WRITE(G, 7 VANGLE (3) 
7 ORNATE Owe De ieee ear oss) 
AB) sh (2) fSUNVD23) 
WRITE(6,8)A(3) 
Bee ROR MA oO nea hies at eeo eh 27°} 
DD(2)=D23+DIP(3) 
UU(2)=D23-DIP(3) 
ZOD AEC OS NO )--COS(UU l2)9) }) 
Seip cau Nae CCOS( DDI 1) )+COS(UUCTI I ))) 
WRITE(6,9)ZD(2) 
OPPORN AI Ole 70 eeee.8, KM" / ) 
ZO ONOH COs (Op eo i=COos( UU 2) )))) 
eee tii ZU ie COS DD 1) JF COSTUUI1))))} 


Wave, ie Ty as 


ons a10au ava 32) 
eat} 130 )3V THs 
-FIATS AIBH 3 
, et 
; (orvera, (orate Dr SUIA, 
" POL OS Ai Ie aby: 
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WRI EA Of 1 Oo eZ UN OF 
PORMAT Oeer® ZU Fae | POetGs TKN oe 


CePGENER AE 7A ONS Ona say PS 


E+ARSIN(A 


DO 100 N=4,NTOTAL 
Tak = 1) DT Paper ARS TN CARETS) 
U(1)=DIP( (2) +ARSIN(A (1)/VU(N 
INDEX=N-2 
DO 200 J=2, INDEX 
DL(J-1)=ARSIN(A A(u 
Be (ai EAR SONA 
J 
TP 


/VD(N)) 
) ) 


DU) =DIP Cu) — hace 
UE cy es) an kd co i Mag 
CONTINUE 
D(N-1)=0. S*(ARSIN(AL 
earn 22) 


Mp 
)/ 
aa 
(J 


N- 1 I 

(N *SI ) | 

DIeEN = DIP(N -1)+0.5*(AR /A(N-2)*SIN(D(N-2))) 
E-ARSIN(A(N Wining 2)*SI ) 

ANGLE (N) =57. 2958*DIP(N) 

A(N)=A(N-1)/SIN(D(N-1)) 
WRITE(6,11)ANGLE(N),A(N 

EORNA EGLO bear OS 

SUMD=0.0 

SUMU=0 .-0 

DO 300 JJ=1,INDEX 
DD( uJ) = Tear 
UU (ud) =U(dd)-DIP 
sat Renee ) 
SUMU=SUMU+ZU (Ju) 
CONTINUE 


PRCMO PE VELGCINY=’ -Fo. 2, KM/SEC’ 77) 


~~ MN ne 


CONTINUE 
Sor 
END 
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19 
G THEF RAW TRACING PROGRAM. 


DIMENSION THETA(46),H(46),V(46) ,PH1(46),PPHI(46) ,UP( 46) 
DEMENSTON TOABEIR( 10).,"R FDS MeS0G)} IP RESUMI(.3:00 } , Peat Clee ee 2 
DEMENSTON TRO IST (300) SRR TSUMC3000), 7TTSUM 0300), 46 ) 

Co er DSi, RRa SUN SUN: ARE. THE FINAL OUGPUT ARRAYS 
DIMENSION CHI1(46) ,CCHI ( pee BEE Lea SR mr due 
DIMENSION D(46), CC 46), D(46), X(46),Y(46),XX(46),YY(46) 
DIMENSION Evy Gi, PATH( 46). TPATH(46) , PPATH( 46) , TDRATH( 4 6) 


NSURF IS THE NUMBER OF INTERFACES IN THE MODEL INCLUDING 

THE FREE SURFACE.OFFSET IS THE SHOT TO SHOT DISTANCE. 

THE CONVENIENT SYSTEM OF COORDINATES IS:POSITIVE X AXIS 
POINTING DOWNWARDS AND POSITIVE Y AXIS POINTING AT THE 
SHOT.THE ORIGIN IS AWAY FROM THE SHOT AT THE OTHER END 
OF THE PROFILE. . 

IF THE SHALLOW PART OF AN INTERFACE IS UNDER THE SHOT 
THE DIP ANGLE OF THE INTERFACE IS POSITIVE, OTHERWHISE 
IT IS NEGATIVE. 
SHOT=1.0 FOR THE RIGHT SIDE SHOT OF THE PROFILE.SHOT=-1.0 
FOR THE LEFT SHOT. 
41 FORMAT (S9A4) 
READ(5,12)TITLE 
12 FORMAT(2A4) 
READ(5,3)NSURF ,OFFSET,RVEL, SHOT 
READ(5,47)LABEL 
47 FORMAT ( 10A4) 
READ(5,41)LABEL1 
C ONE INCH IS TEN KM IN THE RAY PLOT 
SCALE=10.0 
C LEFT SIDE DEPTHS,HH USED ONLY FOR PLOTING THE MODEL 
SULEORMAT CLSSh7723) 
READ(5,331) (THETA 
READ(5,331)(V(1 
READ(5,331) (UP( 
READ(5,331)(H(1 
READ(5, 331) (HH( 
331 FORMAT(17F5.2) 
C CONVERT THE DIP ANGL 
DODROMEN X41 24,7 
THETA(NX)=THETA(NX)*0.017453 
2 CONTINUE 
CALL PLOTS 
C PLOT THE MODEL 
CALL PLOT(2 
CA ULEEPIOI-O 
CALL PLOT(O 
CALL PRLOE(O 
CALL PLOT(O 
CAI RIO UOe 
fj 
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EB, -om0 724) 
GALIIS MBO A CRS Eades 
CALL SYMBO 
GALL Ager Se 
GAWLTAES 2 
CALA PLO 
Cl LaiNre rs 
DO 600 IB= 


QO. 0 

de 4, TARTEES Oe 
DEPTH IN KM! , ; 0. 0, Orns a 
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AILEW 05105 Sy) 
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b/c Be. 0, (Oh rei Rehs 


,9.t-=TOHe. 2409GRR BHT) Ae ORS aC at's 
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(He YY aaie 
(LSRIHTASST Cad 
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JOMAT?S £0 
IMA. A AVLiT CEGRtS Be tees 
shi TA BALTMIOS. 21k + DyEtte Ry . 
Ov3 A3HTO SH? TA: Pon mt wane AW: 
(UR , a +.) 
Jc [HWA SHI ; Vi 


Ae 
am ti i 


42, /BVa, THON, a, uh An 
7 VIGAIE 
so Ate 
rJaRAs rh 
7049 VEX GAY AEWA ATE et 


gt} 
2)9 
3) th, 
AM 


Tenelee Ee Cmeer 10 G00.” 130 
Peerieoe ice TO 600 
Chi Ree ESET / SCALE 0.0.2) 
Celie CU nGerSET/ SCALE “NUP GieN 
Cee ce Mensa rit [Bay SCALE 2 
CAL leIPLO euonUreteac. 3) 

600 CONTINUE 


et 


GAL Dan On uOnms bale) SCALED.O0l27 
CALs iO Te Opis el/ SCALES] (UP G12) ;SCALE), 3) 
CALL, PLOTUOE OS ee Act DAASGAGE 2) 
CAs PaO. 0.0.0; 3) 
CALLS PLOte OF ESET SCALE 0.072) 
CALE. PE OITAGr BCR aCae re GU Psi a BAS CALE) 633) 
CALLS PL OA 0: Os a Hite 7)/SCALE), 2) 
CARI PLO MsOs Os.002.0' 639) 
EMH Os ak 20", CALL Pegi crore SCALE 0.0.=3) 
C WPLOledHe, LRAMECIRLME «FRAME 
CAL Ls PLOW 0p5%9240 jee 
CAD wAXLSA OmOe 00° DISTANCE IN KM! Pee wer Orme Us, Or 
E0205 37, 5) 
C FOR THE DOWNDIP SOURCE MAKE THE T AXIS FOUR INCHES LONGER 


Chale eno OnUr ine NeSee! 2 it. 805,.0.070".0,.1.5) 
CALMS Si MBOL Glen 0p-b., Onei4 ABE ly, 0:10;;40:) 
Ni=4 
N2=NSURF 
301 CONTINUE 
DO S00 N=N1,N2 
[FAN ~EQ26) GOs d 0, 900 
IFAN. EQ: 71 GOg4 0) 900 


COMEACIORSCONTERG@ESeIHE DENSITY OF THE TRACING. FOR N=4 
C YOU NEED HIGH DENSITY.FOR THE OTHER INTERFACES THE 
Ce DENST DY IS yEQWERED: 
FACTOR=0.20 
IF(N.NE.4)FACTOR=0.3 
TEUNCEQ G74) FACTOR=0. 4 
WRITE (6,44)N,LABEL 
44” FORMAL (i eis PEND EREACES 455 OA4/) 
AS uRORM Cd far. 4) ORESE T REDTIME TIME’ /) 
AQe RORMATA © eiien KM SEG SACcUN 
X(1)=0.0 
Y(1) cOEF SET 
A(1)=0.0 
B(1)=0.0 
C(1)=0.0 
D(1)=0.0 
C DEFINE sIBESEQUATIONGP@raelhE, INTERFACES. 
GaOTICE pion) eeVEN hee INTERFACE 1S NOT PLOTTED IN THE 
C RAY DIAGRAM,OR- THE TRAVEL TIME PLOT, THE CORRESPONDING LAYER 
C IS TAKEN SENT OSCRORERSAGCOUNT IN THE CALCULATIONS. 


DO=10C0 Ml =2eN 
Atha, TAN ee aor 4 
B(1)=OFFSET+(H(1)/TAN(THETA(1))) 
1000 CONTINUE 
C CALCULATE THE MINIMUM TAKEOFF ANGLE ZETA(1) MEASURED WITH 
C RESPECT TO THE VERTICAL TO THE FREE SURFACE. 
C THE CALCULATION IS NEEDED IN ORDER TO AVOID ANY UNDEFINED 
C ANGLES ALONG THE PATH OF THE RAY. 
Dec Sel hoe. LO 
RATIO=V(N-IVXX)/V(N-IVXX+1) 
WRITE(6,778)V(N-IVXX),V(N-IVXX4+1),RATIO 
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SW OUTAJUOJAD 3HT WL PRMGSOA SEaORGD Al M3: 
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FORMAT(3F6.2) 131 
Lo tper ica. |) GOTO sts 

CONTINUE 

CONTINUE 

ANGLE=ARSIN(RATIO) 

WRITE(6, 688) IVXX 

FORMAT(15) 

DO 333 IAEK=1,50 

IF(N.NE.8) GO TO 894 

ANGLE =ANGLE+ABS( THETA(N-IVXX+1-1AEK) )-ABS( THETA(N-IVXX 


D+2-IAEK) ) 


GO TO 895 

ANGLE=ANGLE+( THETA(N-IVXX+1-IAEK) )-(THETA(N-IVXX+2-IAEK) ) 

ANGLE=ARSIN(SIN( ANGLE) *V(N-IVXX-IAEK)/V(N-IVXX+1-I1AEK) ) 
F((N-IVXX-IAEK).EQ.1)GO TO 278 

CONTINUE 

ZeT i 2909 1o=000) 4254 ANGLE+THETA(2) 

DES) (=ZEml 154.296 

WRITE(6,8)TEST1 


( 
ZETA(1) 2ZETAL 1) -(FACTOR*0.01745) 
WRITE(6, 45) 
WRITE(6, 49) 
DORN (ieewe meee eke how, Gedy. DEGREES 4 /)) 
IR=0 


DOy 100m te mee 5 

Leuba pheZe Po at Ge OGLOR*O yO)? 45) 
TEST=ZETA(1)*57.296 

WRITE(6,26) TEST 
IF(TEST.GE.90.)GO TO 999 

FORMAT (F7.3) 


Caper ING sinew et RS eOge ker Lec hLON AND REFRACTION ANGLES 


C 


COO GC) © 


DiQi@i@@ 


(PHI,CHI),(PPHI,CCHI) BY USE OF SNELL’S LAW. 
PH1(1)=90.*0.01745-ZETA(1)+THETA(2) 
CH1(1)=ARSIN(SIN(PHI(1))*V(2)/V(1)) 
INDEX1=N-1 

DOS 20eK=2 -INDERI 

PHI (K)=CHI(K-1)+THETA(K+1)-THETA(K) 


AVOIDSCHIAN  tjeWalChe LS NOT DERINED FOR OVERCRITLCAL 
Poa New sole bee dal) plomCALCULATEDTONLY ‘PHI(N=1) 


tS 


ALBOWED: TOBE OVERCRITICAL.. 


THERE IS NO OVERCRITICAL ANGLE FOR A RAY INDIDENT AT THE 
TOPRORMASEVZ 


220 


S20 


IF(K.EQ.INDEX1)GO TO 220 
CHI(K)=ARSIN(SIN(PHI(K) )*V(K+1)/V(K) ) 
CONTINUE 

PPHI ( INDEX1)=PHI(INDEX1) 

L=N 

INDEX2=N-2 

DO 320 KK=1, INDEX2 

L=L-1 

CCHI(L-1)=PPHI(L)+THETA(L+1)-THETA( L) | 
PPHI(L-1)=ARSIN(SIN(CCHI(L-1))*V(L-1)/V(L)) 
CONTINUE 


DEF Mees ee COUATIONS OF TRE INCIDENT AND»REFLECTED RAYS 


C 


eer ee oy whee NCHOEN] »@C.0D FOR THE REFLECTED RAYS. 


RS teceeCUbAL Ea HEE COORDINATES OF THE POINTS AT WHICH 
Tee OUP EN i Nb erePLeCIeD RAYS INTERSECT ‘THE INTERFACES 
Byam yt NGe THe SioueMeOR [HE-LINEAR EQUATIONS. 


Chetan o = (TPT ANCZETALT))) 
D(2)=OFFSET 


7e4) y ; 7 
7 xvi bay - J { Ti28 
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Y. i 


( (NFA -S+ wi-wh 
({MSAT=<24xe 


sarin rat a 


a¢T TA TM30LOML Yaa & spy one 


@) GDie 
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SZ 


ise, 


)+THETA(LM-1))) 
-2)-THETA(LM-1))) 


)/(C(LM)-A(LM) } 
420 CONT 


PPHI(LL)+THETA(LL+1)) 
L+1)-(XX(LL4+1)*TAN(PPHI(LL)+THETA(LL+1))) 
12 2A ae Ed Pay yO Ee 
*DD St ey ee cl Ea ery a Bake 
50 CON 
Ys 1) 
any 1S ae ones OF THE RAY WITH TAKEOFF ANGLE ZETA(1) 
THE LAST REFLECTING HORIZON OF THE RAY IS INTERFACE N 
FDIST HAS TO BE AN ARRAY .ALSO RRTSUM,TTSUM 
FYY(1)=OFFSET-YY(1) 
INCREASE eZETA CI) 4IG@aDECREASE THE OFFSET 
IF(ABS(FYY(1)).GT.OFFSET)GO TO 100 
TP CBG lei bsOs)GOet0 999 
IR=IR+1 
EDT SieGCTR = ey) 
CALCULATE THE TRAVEL TIME OF THE RAY 
INDEX1=N- 1 
DO 60 NA=1, INDEX14 
PATH(NA)=SQRT((X(NA+1)- 
TPATH(NA)=PATH(NA)/V(NA 
60 CONTINUE 
DO 70 NN=1, INDEX1 
PPATH(NN) =SQRT ( (XX ( 
TPPATH(NN) =PPATH(NN 
70 CONTINUE 
TSUM=0.0 
DO 80 IP=1,INDEX1 
TSUM=TSUM+TPATH(1IP)+TPPATH(IP) 
TTSUM(IR)=TSUM 
80 CONTINUE 
RTSUM=TSUM- (FYY(1)/RVEL) 
RRTSUM(IR)=RTSUM 
CONS T=18-0£0..04745 


AE A) )**2+(Y(NA+1)-Y (NA) )**2) 


NN+ 1 
)/V 4 


N Ge rN ers Oe aie 
/V 


NN 


PRO OU Me bne eke mee i ob eOr Tee HORTZON N 
DEAN CNS ar eG VEN ome, 10-18 
IF(ABS(PHI(N-1)-ARSIN(V(N-1)/V(N))).LE.CONST)GO TO 52 
(SO WRITE ODO YY tink SUM, TSUM 
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CALL PLOT(SHOT*OFFSET/SCALE,0.0,3) 

BONG Ot aac N 

CALE PLOmSHOTRY (CIAIYSCAUE =X CIA) /SCALE 2) 
CONTINUE : 

LX=N 

DOL602eetD=diENDEX4 

LS 

GRU SPIOmNSRO Tay lek SCALE Y= XL) SCALE, 2) 
CONTINUE 

CONTINUE 

CONTINUE 

WRITE G27 }IR 

DO 800 IL=1,I1R 

FFDIST(IL)=FDIST(IR-1L+1) 

CONTINUE 

DOmS Oar Bean eR 

FRTSUM(ILL)=RRTSUM(IR-ILL+1) 

CONTINUE 
FFDIST(IR+1)=0 
FFOL Sih Re) = 
FRTSUM(IR+1)= 
FRTSUM(IR+2) =1 
CALL LINE(FRT 
CONTINUE 

CALLE PL OTeO 208 Ov 04999) 
STOP 
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